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ABSTRACT 


Numerical models for nearshore circulation patterns in the surf zone 
have been developed and applied to an observed condition subjected to 
a sea breeze environment. Bottom topography and input waves were 
derived from observed data to predict surf zone circulation as a 
function of time of day. It was found that many features observed in 
the surf zone were modeled but wave-current interactions are known 


to be important. 


Wave-current interactions were modeled for shallow water assuming 


a two-dimensional motion which included rip current and longshore 
current components, The refraction effects caused by even small 
currents produce major changes in the wave induced driving forces 
in the surf zone which leads to the prediction of entirely different 
rip-current patterns when wave-current interactions are considered. 
Numerical results are presented and a discussior ut the numerical 


techniques is included. 


A review of water wave theories t> include mass transport, vorticity 
and current was made for a vertical section in shallow water of 


constant depth. 


—— 
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l. INTRODUCTION AND SUMMARY : 
- 4 bel INTRODUCTION 


. In the nearshore area waves arriving from offshcxe continuously 


bring in momentum, energy and mass. Since the shoreline provides 


a fixed boundary the momentum and energy fluxes are dissipated in 
the surf zone. Most of the energy is converted to turbulence in the 
breaker zone but enough is left to supply a nearshore current 
system and move locse bed material. The momentum brought in by 
vet 4 the waves will drive the nearshore current system and cause a local 


set-up or set-down of the mean water level. 


+ Over the past four and one half years a series of analytic developments 
have been attempted to model some of the more pertinent character- 
istics of the surf zone. The work is conveniently divided into two 
: broad groups: statistical and deterministic. In the statistical 
oe approach (Collins, 1971 and Collins and Wier, 1969) a relatively 
simple beach topography was assumed and the effects on wave 
height statistics computed together with longshore currents and 
wave set-up. More recently, (Noda, 1972, 1973) a deterministic 
approach employing monochromatic waves and much more cor. plex 


is beach topographies has been explored. 


4 number of sub-tasks have been investigated during the past year. 
The three specific sub-tasks receiving intensive investigation 
include: 
a) the application of wave-induced circulation computations 
on beaches having rythmic topography. 
b) the development of a numerical model for wave 
induced circulation which includes wave-current 
interaction. 


c) the analytical investigation of wave, current, and 


vorticity interaction. 


PRR DRE tee oo 


The following sections of this report present details of the work 
performed. The subsections below pre 2nt a brief review of some 
earlier work and a summary of the work completed during the 


past year. 
1.2 REVIEW OF EARLIER WORK 


In a recent study by Noda (1972, 1973) the solution to wave-induced 
nearshore circulation due to the incoming wave-bottom topography 
interaction vas studied. Results for both normal and oblique wave 
incidence were presented and while the results generally agreed 
with recent field data from Sonu (1972), Figures 1.1, 1.2 and 1.3, 
the numerically derived circulation velocities tended to be larger 


than measured in the field, Figures 1.4, 1.5 and 1.6. 


Several possible reasons for the apparent discrepancies can be 


postulated including; 


a) neglect of wave-current interaction 
) bottom fiction approximation 
c) chuice of wave breaking criteria 
d) assumption of monochromatic waves which conse- 


quently all break at the same location 
e) over-estimates of the incoming wave height or 
errors ic. direction 


f) approximav.ons made in the analytical developments. 


Of the possible reasons for dif“cr-uces the assumptions made to 

comply with (c), (d) and fe) srcduce simiiar effects in that the 
nearshore circulation pat.ern is strongiy influenced by the wave 
breaker location. The dominant driving forces are produced by the 
radiation stresses induced by breaking waves. Also, because of 

this it must be realized thit «en relatively weak currents change 

the breaker location and characteristics hence, the importance of wave- 
current interaction ig 1 «.ajor one. Therefore reason (a) is of prime 


importance, 
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3. LONGSHORE SURF ZONE PROFILE 


Figure 1.1. Dependence of Current Patterns on Wave Incidence 
j Angles and Surf Zone Topography [From Sonu, )972]} 
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it is believed that bottom friction effects are important and hence 
are never ignored in these investigations. However, approximations 
are necessary in order to yield a tractible numerical model. I is 
apparent that considerable room for improvement exists in the 
approximations generally made. 


The numerical techniques employed were found to influence the 
predicted currents and a certain effort was needed to refine the 
earlier more crude methods. New approximations include the choice 
of rnore realistic beach topographic model and the procedures to 


solve the governing differential equaticns. 


The following subsection (1.3) of this report presents a brief 
summary of the technical work which has been oriented towards 
improvements and refinements in the nearshore circulation models, 
More complete details are presented in Sections 2 and 3. 


1.3 SUMMARY 


1.3.1 Wave Induced Circulation Over a Rhythmic Topography 


The data obtained bythe Coastal Studies Institute (SALIS by Sonu 

et. al., 1973) has been investigated and attempts have been made 

to model the wave induced circulation using the procedures developed 
by Noda (1972) in an earlier phase of the work. 


The steps required are: 
a) topographic model 


b) wave height-wave direction field 
c) solution of the momentum equations 
d) comparison with observed data 


The field data was used to provide (a) and the offshore wave conditions. 
The topographic model was developed by choosing empirical functions 
and constants to closely simulate the observed topography. Fig. 1.7 
presents 2 sample of the topographic simulation as used. More 
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Figure 1.7: Comparison of Analytical Model and an Observed 
Nearshore Topography 
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details will be included in Section 2, 


The wave height-direction field has been computed using the ray 

equations (see Noda, 1972) but some variations on the numerical 

techniques were found which yielded significant improvements in 

accuracy and speed of computation. The revised technique is 

based on a relaxation procedure rather than the previous method of 

marching along rays. Wave heights and directicns are computed 

directly at the required grid points and yield considerable savings in 
computing time otherwise spent on interpolation subroutines. The techniques 
are fully detailed in Section 3. 


The solution of the momentum equations for the wave induced circu- 
lation follows the pruccdures outlined in Technical Report No. 3 
(Noda, 1972). 


One apparent difference between the numerical results and typical 
observed data is that the numerical predictions show a too strong 
concentrating effect of the circulation near the transverse bar 

(or shoal) and a number of localized eddies in the nearshore area. 
Possible reasons for such effects have been indicated in Section 1. 2 
and these are alsu discussed in Section 2.4. However, in spite of 
some obvious shortcomings it is apparent that the approach and 
results outlined in Section 2 have yielded a reasonable modeling 
capability for many features of the nearshore circulation over a 
rhythmic topography. 


1.3.2 Wave and Current Interaction Over a Rhythmic Topography 


It is apparent from even a casual glance at a beach that incomitg 
waves interact strongly with the local currents which are themselves 
induced by the waves. Section 3 presents a detailed analysis of 

this problem. The interaction produces two dominant effects; 


li 


wet csualil 


the currents change the wave refraction and also change the breaker 


locations, 


These wave-induced nearshore circulation patterns were derived 
assuming no wave-current interaction. Thus interest was developed 
to determine if the effects of wave-current intc “action produced 
signi‘icant changes in the nearshore circulation 2atterns as observed 
in prototype. Section 3 deals with the theoretical development and 
numerical computations of this process as affecting the circulation 


within the nearshore zone. 


The initial computational steps followed those given in Section 1.3.1, 
i,e. topography-wave height and direction field-solution of the 
momentum equations. Thenthe resulting circulation ‘relocities 
were considered as an existing mean current system, and waves were 
again propagated into this system and a new wave height, direction 
and nearshore circulation pattern obtained. It was hoped that con- 
tinual interaction would finally yield an "equilibrium" solution in- 
cluding wave-current interaction, However, attempts to directly 
impose this derived mean current system in an interaction process 
with the incoming waves lead to failures of the technique because the 
mean current system derived for no wave-current interaction was 
too large. Hence, conditions arise where the local waves were no 


longer able to propagate into some areas, 


An attempt was 1nade to take only a percentage of the initially 
derived current system and then inc!vde interaction with the waves. This was 
partly successful and indicated that considerations of wave-current 
interaction were extremely important. Some major changes in the 


computed nearshore circulation systern were produced. Section 


ee ee eT 


3.3.3 presents some of the results. 
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it has been demonstrated that wave-current interactions are of mzjor 
importance in deterrnining the nearshore circulation but a complete 
solution was not possible because of the occurrance of regiona in the 
nearshore zone where waves could no longer propogate when oppo sed 
by a current. Two major conclusion are deduced, 


a) the current-wave interaction theory needs further 
development to include the special case of ''no wave 
propagation" in some regions; 

b) the nearshore circulation system is basically a 
non-steady pulsating system in that the breaking 
waves initially produce a circulation system which 
shuts off the waves in some regions and decays 
until the waves are re-established and reproduce the 


initial circulation. 


There seems to be a considerable amount of qualitative field data 


to support the second hypothesis (see, for instance, Sonu, 1972). 


1.3.3 Wave-Current Interaction with Vorticity (Two-Dimensional) 


As waves approach a shoreline they transport energy, momentum 
and mass from deep water towards the shore. Many aspects of 
the momentum balance have been evaluated in Sections 2 and 3 and 
summarized above. Mass transport by waves is closely related to 
the vorticity present in the water column. Section 4 of this report 
presents a detailed review of wave motion, currents, mass transport 
and vorticity and their interaction in two-dimensions. The early 
work of Dubreil-Jacotin (1934) and others is reviewed. It is shown 
that there are an infinite number of solutions for periodic waves in 
a perfect inviscid fluid associated with the presence of a more or 
less arbitrary vorticity distribution. A current having a velocity 
profile which varies over a vertical has an associated vorticity 
distribution and hence the form of periodic waves present do not 
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necessarily follow the classical Stokes solution, 


In Section 4 of this report the equations required to solve at least up 

to the first order the problem of small amplitude weve propagation in 
the presence of an arbitrary current for an arbitrary wave spectrum 
have been presented, and a review of the special solutions previously 
obtained for a single wave length has been made. The problem in 


general requires lengthy numerical computations. 


However, for a current whose velocity distribution can be approx- 
imated by a linear depth dependence, it has been shown that, at 
most, a sincle numerical quadrature was required to obtain the 
average velocity components. This method may then be used to 
estimate the forces due to wave action in the presence of a current. 
An experimental knowledge of the current velocity at but a few depths 
(two minimum) will define the parameters neccessary to completely 


solve this problem, 


2. CIRCULATIONS UNDER THE SEA BREEZE CONDITION 
2.1 INTRODUCTION 


When the nearshore wave field is strongly influenced by a sea breeze, 
local wind waves undergo diurnal changes in height, period, and 
incidence angles. In the northern hemisphere, the wave direction 
rotates clockwise, while heights and periods both grow steadily 
toward late afternoon. Usually, a background swell is superimposed 


on these wind waves. 


Nearshore circulations, which are sensitive to breakers and their 
incidence angles, will undergo rapid changes accordingly. Diurnal 
changes in nearshore and surf zone topography under this condition 
are probably more gentle. This situation is known to develop at a 


number of tropical and subtropical regiong of the world. 


In this chapter, a series of computations are performed to simulate 
successive stages of nearshore circulation under the influence of a 
day-time sea breeze condition. Some of the basic considerations 


included in the present computation are surnmarized as follows: 


1) In reality, the change in the circulation velocity field occurs 


wales ETE weidee oS 


as a continuous process. However, a finite-difference solution of 
time-dependent equations involves technical difficulties as well as a i 
considerable amount of computer time. Instead, the computation is 


performed for four discreive stages of circulation development (at 


three hourly intervals) using steady-state equations. 


Vint q 2) Quadratic inertia terms impose difficult, if not insurmount- 
able, restrictions to the computation. Consequently, the equations 


of motion are linearized by neglecting the inertia terms. 


3) Velgcity variations over a vertical are neglected. 
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4) The foriaulation of the bottom friction term in the momen- i 
tum equations was derived following the assumption that circulation 
velocity components are small as cornpared to wave orbital velocity, 
as in the previous report (Noda, 1972; Thornton, 1969). 


5) The effect of interactions between wave and circulation, 
as discussed in detail in Section 3, is not included in the computa- 


tions presented in this section. 


6) When wind waves and swell coexist as separate wave trains, 
there will be an interaction not only between them but also between 
the currents they drive simultaneously. This situation is extremely 
complex and involves a number of mechanisms wnich are not well 
understood. As an alternative, the case of coexisting wind wave and 
swell is treated by vector addition of the velocity fields associated 


with each of the wave trains. 


2.2 GOVERNING EQUATIONS 


The method of computation is to solve by a finite difference approxi- 

mation a set of steady-state linear equations of motion and a continu- i 
ity equation. Basic mathematics of this method have been discussed 

in detail in the previous report (Noda, 1972). However, for the 


benefit of the reader, these will be briefly summarized: 


Equations of mction (vertically integrated) are: 


an _ = , 
8 3x ae a; ees 
SNe a 2 2.2 


8 ay y y 


and a continuity equation is: 


& Lunta)]+ 2 [vinta] o. 2.3 


where x and y are taken normal and parallel to the coast, res- 
pectively. 


M, and My denote radiation stress terms (Longuet-Higgins, 
1964), given by 


, 1 Bx | OT ey 2.4 
Ne ~~ 9 (n+d) ax t ay ) 

2 Na ae Se 
My - pinta) Vay * Ox ) 2.5 


where, in shallow water, 


ar 2 [ 2 er ] 
Ox te 98H 3 cos @+ sin’ 6 2.6 
sail 2 [ .2 2.) 
Cy ie Se 3 sin 6 + ccs 8} 2.7 
! 
and ; 
_ ol 2422 
Teg ge Te pgH sin 6 2.8 
The friction terms are simplified as: 
- 2c Hu = 
F. *(qFd)T sinhkd * F.d.u 2.9 : 
F cou = Fidv 2.10 


y  (n+d)T sinh kd 


where cis friction coefficient (0.01 in our computation); d is the 
water depth, and 1 is a set-up or set-down relative to the mean sea 


level. 
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Defining a stream function given by 


3 = -ud, ot = +vd 2.11 
and assuming 
nide=d 2. 12 


Equations 2. 1-2. 3 reduce to a single equation: 


oF oF 


avy , ay , dy 8a, ox By 
2 2 F 3éay F ax 


a.) 3 - 3 3 
4 {2 (2 Gee @))- 2G S| 0 a 


The boundary conditions are: 


Se = 0 at x=0 and @ , 2.14 
and 
tly,z)=¥ (yt), x) 2.15 


The latter condition implies that the circulation field is periodic 
along the shore at a spacing equal to the wavelength 1 of the 
bo'tom topography. 


The computation solves Eq. 2,13 using a relaxation (or Gauss - 
Seidell) method, as already discussed in the previous report. The 

W values at the inshore and offshore boundaries can be chosen 
arbitrarily. I. this case, is chosen to be zero at x=Oande. 
The iterative procedure was continued until a condition 


ly. 


j+1 7 9,/714,1 s 0.05 2. 14 
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was achieved between successive interation cycles LF and *. 41° 


The radiation stress field to be entered into Equation 2.13 is provid- 


ed from a wave ray equation which combines effects of shoaling 


and refraction 


2 
PE + pls) PE + alsin- 0 2.15 
Ds . 


where 
p(s) = -cos of 22] Beer 2 ac] 
2.16 
2 2 2 
4 1 Cc 2 1 Cc 
q(s) = win’ E xe seine cos ole ac] +cos de ae] 
ax ay 
where 


s is the arc length along the ray 


B is the wave intensity, and C is the celerity. 


Previously, the ray equation was solved by a fourth order Runge- 
Kutta scheme. In the present report, this equation is sclved by a 
relaxation technique, as described in detail in Section 3. This 
method computes incident w..ve heights and angles directly on the 
grid, whereas the previous method traced wave rays individually, 
which required additional visual inspection of wave ray density and 
interpolation steps to transfer the ray data onto the grid. The new 
method thus allows the entire wave field computation to be carried 
out in a single run of computer processing, resulting in a substan- 


tial improvement with respect to both speed and accuracy. 
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2.3 BOTTOM TOPOGRAPHY 


The input information for bottom topography and wave characteristics 
is derived from the observations carried out by the Coastal Studies 
Institute on Santa Kosa Island, Florida, in 1972 (SALIS Project, see 
Sonu et al., 1973). The CSI data are especially pertinent ic our study 
because they contained detailed characteristics of the surf zone tops- 
graphy to which the nearshore circulation is known to be sensitive 
(Sonu, 1972, 1973). The CSI data aleo contained general information 
cf circulation pattern and current velocities as revealed from re- 


peated dye experiments. 


Importance of bottom topography, particularly that of undulations in 
the surf zone bottom, to nearshore circulation has been pointed out 
by a number of field observers, among them Evans (1939}, McKenzie 7 
(1958), Shadrin (1961), Davis and Fox (1971, 1972), and Sonu (1972, 
1973), Surf zone topographies as reported by these investigators E 
are summarized in Figure 2.1. Evans reported a meandering current 1 
consisting of an inflow acroas the bar and an outflow originating from , 
the shoreline embayment, McKenzie reported an inflow across the : 
bank (shoal) and an outflow along a conspicuous rip channel (or de- 
pression) between banks. It should be noted that although a schematic 
presented by McKenzie depicts the shoreline with a straight line, his 
photographs indicated a periodically curved shoreline. According to 
Shadrin, an outflow generally initiated in the embayment, but its 
orientation depended upon not only wave direction but also wave — 
height. Davis and Fox reported meandering currents under wind 
wave conditions, These rhythmic topographies had wavelengths rang- 


ing between 70 and 200 meters. 


Figure 2.2 shows the surf zone topography at the site of the CSI pro- 
ject. Note that a cuspate portion of the rhythmic shoreline descends 
directly to a shoal in the surf zone. A line of iongshore bar exists 


approximately 30 meters from the average shoreline position, This 
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rhythmic topography was formed at the time of a strong local storm 
and remained essentially unchanged for as long as 14 days while a 


local sea breeze dominated the area, 


Figure 2.2 also shows a typical example of water movement as re- 
vealed from the movement of dye. The dye, initially injected at the 
break point on a shoal, streaked toward an embayment in approximately 
the same direction as the breaker. It then travelled parallel to the 
shoreline for some distance before making a seawardturn. The out- 
flow across the surf zone usually occurred on the depression, which, 
upon reaching a break point, tended to turn alongshore and eventually 
returned shoreward across the downstream shoal. This type of mean- 
dering current pattern was typical of afternoon conditions when wind 
waves associated with the sea breeze arrived obliquely to the coast. 
Current speeds in the meandering currents generally amounted to 

30 cm/sec in the inflow current across the shoal, 10 - 15 cm/sec in 
the parallel current near the shoreline, and about 20 cm/sec in the 


outflow or rip. 


During the morning hours when the wave field was dominated by the 
background swell, the currents tended to form closed circulations 
of minor velocities, consisting of an inflow on the shoal and an out- 
flow on the depression, Maximum speed under this condition was no 


more than 20 cm/sec, 


For mathematical representation, a rhythmic topography may be 
broken up into three components, (1) mean profile, (2) longshore 


bar, and (3) longshore undulations, 
The mean profile of a c”~st is generally concave upward and may be 
approximated by 


qd, = ux’ 2.17 
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in which qd) is the depth measured from the mean sea level, x ig 
the distance seaward from the shoreline, and y and y are numer- 
4 ical coefficients; especially, y <1 to ensure the concavity of the 
_ 4 : profile, Bruun (1973) shower on the basis of a wide range of evidence, 
that y varies between about 2/3 nearshore and about 1/2 offshore, 


The bar can be defined, for the sake of simplicity, as a symmetrical 
hump superimposed on the mean profile. Assuming a bell-shaped 


configuration similar to an error function, the bar profile is given 


by, 


d, =b. exp [- “x, )?H(x,/2?"| 2.18 


The longshore undulation is generally confined within the surf zone, 
and its amplitude attenuates quite rapidly outside the breaker line. 
Thus, we assume a longshore undulation whose amplitude decreases 


linearly toward zero at x = ly i.e. 
d, =a (1-x/1,)sin 2 (y-6) 2.19 
in which a is the maximum amplitude and \ is the wave length of 


the undulation, The term § in Equation 2.19 represents a degree 
of distortion to be introduced in the geometry of the undulation, 
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Normally, this will consist of two parts: 


b= 6, + 6. 2.20 
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Where a longshore current is significant, the longshore cross-section 
of the sinusoidal undulation is skewed, yielding a steeper slope facing 
the downstream side. Furthermore, under this condition, the crest~ 
line of the undulation will extend cbliquely seaward from the shore- 


line. 


The first of these effects, the skewness, can be incorporated in 
Equation 2.19 by considering 8) of the form 


. 2n 
8) = binag Sinz (y-8)). 2.21 


In other words, the symmetrical sinusoid of the original undulation, 
sin iy), is distorted by displacing the coordinate y by a variable 
distance §) in such a way as to achieve a steep downstream slope. 
The displacement is maximum (6 ) along the crest of undulation, 
e.g. at y= ut (2nt+1) + 6 , decreasing in both directions away 

pric dé : uae eu) 
from this in proportion tc sin we (y-6,)- 


The oblique downstream orientation of the crest of the undulation can 
be represented by &, of the form 


b5 = xtang 2.22 


in which © is the angle between the normal to the shoreline and the 


crest of undulation. 


Thus, combining the mean profile, a bar, and skewed undulations, 


the general expression for the rhythmic topography is 


d= d,-d,+d, 


= yx’->. exp [-cx-2,)? x, /2”"| +a (1-x/1,) sin = (y-6)- 5) 


2.23 
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Figure 2.3 shows successive superimposition of qd): d, and d,, in 
which «= 0.075, y = 0.600, b= 0.300 (meters), X= 30,00 (meters), 
a = 0, 200 (meters), Ls 80 (meters), } =115 (meters), and a= 20°, 


2.4 WAVES 


Figures 2.4(a), (b), and (c) show diurnal changes in wave character- 
istics. Typically, the waves during the morning were dominated by 
the background swell arriving normal to the shore. As the sea breeze 
began to increase between 1100-}200 hours, small wind waves be- 
came superimposed on swell. Wind waves subsequently grew both 

in height and period, while rotating its direction clockwise, until they 
dominated the sea state around )500-1800 hours in the afternoon. In 
the evening hours after 1800 hours, wind waves steadily attenuated 
and were gradually replaced by the background swell until the next 


morning. 


In Figure 2.4(a) and 2.4(b), it is seen that the wind waves (0.3-0.7 
cps) were strongly coupled with sea breeze, so that the period in-~- 
creased rapidly from about ! sec at 1000 hours to 3 sec at 1600 hours, 
the time of maximum sea breeze. The direction of wind waves also 
increased from about 20° to 40° against the normal to the shoreline 
(Fig. 2.4(c)). The swell spectrum underwent a slight change, its 


direction remaining essentially perpendicular to the shoreline. 


From these data, the wave heights, periods, and directions to be 
input into the computation were determined, as shown in Table 2.1. 


The significant wave height was computed from the power spectrum 


according to 
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TABLE 2.1 


Time of Day Wind Wave 
hours Hiss T 6 


° 
cm sec 


17.9 1.62 0.0 
1500 28.8 2,41 15.0 


33.9 2.75 25.0 


INPUT WAVE CHARACTERISTICS 


26.0 7.00 -2. 
28.8 7.70 -1. 
29.0 7.90 


enn 


} : The term inside the parentheses denotes either wind-wave or swell 
4 portions of the power spectrum, as plotted in Figure 2.4(b), Wave 
: : periods were obtained directly from the spectral density peaks for 


wind-wave and swell. 
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2.5 RESULTS 
2.5.1 Circulations Under Wind Waves 


Figures 2.5 and 2,6 show streamlines caused by wind waves cnly at 
1200, 1500, 1800 and 2100 hours. Note that the streamline separation 


represents .2 m/sec. 


According to Figure 2.5 (case of normal wave incidence), an inflow 

dominates the area of shoals (y= 30-60, 140-170 meters), An inflow 
also occurs at part of the depression immediately to the right of the 

shoal. However, most of the depression area (y= 90-120 meters) is 
dominated by outflow. Thus, there is a general indication that an 


inflow is strong on the shoal and an outflow is strong on the depression, 


However, a detailed streamline distribution is more complex and in- 
cludes some departures from the general rule. There is a small but 
aS well-defined eddy immediately to the right of the shoal (y= 70-90, 
185-205), which surrounds an area marked by a contour 0.4 meters. 
Another eddy of much smaller velocity is located almost directly 
offshore. These eddies have not been noticed during the field -bser- 
vation, It must be noted that, although the congested streamlines 
give the impression of a strong current, they only involve velocities 
on the order of a few cm/sec. Normal velocity components are 
clearly larger than the longshore components. The inflow velocity 
on the shoal is on the order of 1.5 cm/sec; the outflow velocity in 
; the middle of the depression is on the order of 1.8 cm/sec. The 
maximum inflow velocity reaches about 8.6 cm/sec at y= 70; the 
maximum outflow velocity reaches about 7.6 cm/sec at y= 85. 


Maximum velocity outside of the surf zone is only about 4 cm/sec. 


These low velocities are typical of weak breaker activities prior to 
the arrival of the sea breeze wave front in the surf zone. It is 
noticed in Figure 2,5 that waves are breaking only in the immediate 
vicinity of a shoal. 
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Figure 2-5: Streamlines and Bottom Topography for Wind Waves 
at 1200 hours Normal Incidence : 
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Figure 2.6 shows cases of oblique wind waves inthe afternoon. These 
streamlines now exhibit a stronger tendency for meander than in the 
case of normal incidence of Fig. 2.5, The outflow portion of the mean- 
der is located in the depression. However, an inflow also occurs 

at the depression nearer an upstream shoal. Small eddies tend to 
persist throughout the period of computation. 


One of the conspicuous features of the afternoon situations is the 
tendency for the longshore current along the bar crest to intensify 

in proportion to the breaking activity. Fig. 2.7 shows the distribution 
of bveakers at times corresponding to Fig. 2.6, Breaking is the most 
intensive at 1800 hours, e.g. at the peak of sea breeze activity, 
generating a strongest current along the bar crest (maximum 23 
cm/sec). Both before and after this event (e.g., at 1500 and 2100 
hours) when the breaker zone was narrower, current speeds along 
the bar reached a maximum of only 15 cm/sec. Ccncentration of 
longshore current velocity in the breaker zone arises from the 
lungshore wave thrust generated directly by a breaking phenomenon, 
in proportion to the rate of shoreward decrease in the flux of long- 
shore momentum across a plane parallel to the shore. 


2.5.2. Circulations Under Swell 


Figure 2.8 show streamlines associated with swell. Only two cases 
are shown inasmuch as the swell characteristics changed little under 


the sea breeze condition. Streamline separation is .6 m/sec. 


A salient feature of these streamlines is the occurrence of a local : 
circulation immediately to the right of the shoal (y= 60-80, 175-185), . 
which contains velocities ag high ag 120 cm/sec seawards and 80 
cm/sec onshore. These circulation are located somewhat offshore 

of the eddies as noticed in the case of wind waves (compare with 


Figures 2.5 and 2. 6). - 
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Figure 2-6: Streamlines During Oblique Wave Incidences 
in the Afterroon 
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Figure 2-7: Breaker Distribution 
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Figure 2.8: Typical Streamlines Under Swell 
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These velocity patterns contrast strongly with the case of closed 
circulations previously observed by Sonu off the Seagrove beach 
(Sonu, 1972; See Figure 1.2 inthis report). In the latter case, the 
bottom undulation was symmetrical, containing a broad shoal and 
narrcw depression. Inthe present case, a depression occupies a 
larger area than a shoal, and the shoal is non-symrnetrical, causing 
a more complex distribution of radiation stresses than in the case 


of a symmetrical broad shoal, 


2.5.3. Circulations Under Coexisting Wind Waves And Swell 


Figure 2.9 shows superimposition of streamlines associated with 


wind waves and swell. Again, streamline separation is 0.6 m/sec. 


As expected, the results generally indicated both features of wind- 
wave and swell cases. At 1200 hours, when wind waves produce 
weak breakers, the current field is dominated by swell. Effects of 
wind waves steadily increase through 1500 hours toward 1800 hours, 
the tendency for current meander becoming gradually more evident. 
At 1800 hours, a current arriving at a shoal partly escapes seaward 
and partly meanders back shoreward. A local circulation near the 
tip of a shoal persists, reflecting a complicated radiation stress 
distribution over the sharply skewed bottom topography. It is also 
noted that a strong longshore current along the bar crest remains 

in force during the time of maximum sea breeze at 1800 hours. In 
general, current activities are concentrated around the steep fall 


of this shoal where the breaker height variation is most pronounced. 


2.6 DISCUSSIONS 


The simulated streamlines indicate both similarities and differences as 
compared with field observations. In general, the feature of inflow 


dominance over the shoal and outflow dominance over the depression 


is revealed inthe computed streamlines, but it is also disrupted to 


various degrees by the occurrence of localized eddies and srnall 


37 


Figure 2-9: 
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Streamlines Under Combined Effects of 
Wind Waves and Swell 
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circulations persisting near the steep face of the skewed bottom un- 
dulation. Especially in the case of swell, these localized flows 
tend to dominate the overall streamline distribution. Also, the 
computed velocities tend to be higher than observations by a sub- 


stantial margin especially in the case of swell. 


Several approaches seem possible to improve the degree of relia- 


bility of numerical simulation for nearshore circulations. 


Firat, the criterion for breaking inception and the estimation of 


wave heights during breaking should be improved. The present 
computation uses the Miche criterion, 


) = 0.12 tanh 2n (3) 
@ b Ca 


for both breaking inception and post-breaking wave height. Since 
this criterion requires wave heights to diminish to zero at the shore- 
line, the rate of wave height reduction during breaking, hence the 
magnitude of radiation stress, may result in over-estimation. This 


could be one of the causes for overestimation of velocities. 


There exists a critical deficiency of knowledge on the behavior of 
breaking waves. One way to overcome this difficulty may be to 
take into consideration a wave set-up in the water depth estimation 
in the surf zone. This problem has been handled numerically in a 
two-dimensional case (Hwang and Divoky, 1970). Inthe three- 
dimensional case, as in our study, this problem could be handled 
by stepwise approximation, First, the result of the computation 
which is based on the assumption (Eq. 2.12). 


nid=d 


could be substituted into the starting equations 2.1~-2.3 to determine " 
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In the next iteration, n, will be added to the mean-sea-level water 
depth d andthe new wave field and streamlines will be determined, 
This resuit will again be recycled to the starting equations to deter- 
mine nN, and initiate the second iteration, and so on. These procedures 
will result in a slower breaker height reduction on the shoal and 


hence smaller radiation atresses and weaker currents, 


The second approach is to take into consideration the randomness in 
the incident waves. Since wave breaking will occur ina zone instead 
of at a point, the radiation stresses will be spread more broadly, 
resulting in a general lowering of peak current velocities. In the 

case of two-dimensional longshore currents, this approach has result- 
ed in a velocity distribution comparable to a derivation using a tur- 


bulent momentum mixing or eddy viscosity assumption (Collins, 1972). 


Third, a more rigorous formulation of the bottom friction term may 

he needed. Inthe present computation, the bottom friction is asscciat- 
ed primarily with wave orbital motion, Retardation of circulation 
velocity, presumably of considerable magnitude, is not taken into 
consideration in full value. As already mentioned, this approach 
requires readjustment of numerical scheme to ensure a sufficient 


degree of computational stability. 


Fourth, it must be noted that the present computation does not con- 
sider interactions between wave and circulation. Therefore, there 

is an implicit assumption as if the wave field had been abruptly re- 
moved after driving the current instantaneously. However, the current, 
Once produced, will interact with waves at all phases of its develop- 
ment. It is possible that the effect of such interactions is to produce 

an equilibrium circulation with less current velocities than obtained 

in the present computation, or a pulsation of the circulation around 


a certain mean equilibrium state. 
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3. WAVE CURRENT INTERACTION OVER VARIABLE 
TOPOGRAPHY 
3.1 INTRODUCTION AND REVIEW OF HISTORICAL WORK 


The available literature on surface wave-current interaction is not 
extensive. Unna (1942) and Sverdrup (1944) considered the case of 
deep-water waves encountering a following or opposing current and 
applied their results to waves in tidal entrances. Johnson (1941) 
discussed the refraction of deep-water waves encountering a uniform 
current moving at an angle to the wave system, Arthur (1950) studied 
the problem of shallow-water waves being refracted by both changes 
in bottom bathymetry and a nonuniform current system, Application 
of refraction effects due to a current distribution similar to an intense 
rip current was solved by considering the analogous problem of deter- 
mining the minimum flight path of an airplane flying in a variable wind 
field. 


Taylor (1955) investigated the influence an outward flowing surface 
current would have in preventing the passage of waves coming in from 
the sea, This study was in association with the concept of utilizing a 
surface current produced by a curtain of air bubbles as a pneumatic 
breakwater". Evans (1955) performed an experimental investigation 
of this concept. 


Ursell (1960) and Whitham (1960) developed the general geometrical 
equation governing the interaction of a variable current and any type 

of wave motion. Ina classic series of papers by Longuet-Higgins 

and Stewart (1960, 1961, 1962) and by Whitham (1962) the conservation 
equations of mass, momentum and energy per unit area for a wave 
system superimposed on a variable current system were derived. A 
very good summary of this work is given by Phillips (1966). Taylor 
(1962) studied the characteristics of free-standing waves on either a 
contracting or expanding current and provided experimental data. 
Hughes and Stewart (1961) also conducted experimental investigations 


to determine the characteristics of gravity waves on a shear flow. 
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Recently Jonsson, Skougaard and Wang {1970} concentrated attention 
on the ''current-wave set-down" for two-dimensional wave current 
propagation over a gently sloping bed. Kenyon (1971) studied the 
kinemetics of deep-water waves in conjunction with a variable current 
to show the possibility of either the trapping or total reflection of 
waves by the current. 


To carry out the basic objective of this study as indicated at the very 
outset of this introduction, the important kinematic and dynamic 
relationships are first set forth. Numerical techniques are developed 
to solve these relationships so that the stream function and associated 
circulation pattern can be obtained. ‘he basic philosophy is to first 
solve the nearshore wave-induced circulation problem with no wave- 
current interaction. Then the output of these circulation velocities 
are now considered the existing mean current system, and waves are 
again propagated into this system and a new wave height, direction 
and nearshore circulation pattern obtained. It is hoped that this con- 


tinual interaction will finally yield an "equilibrium" solution. 
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3.2 WAVE CURRENT INTERACTION 
3.2.1 Wave Kinematics 


Inherent in the concert of three-dimensional waves is the motion of 
a "wave front’. Crests and troughs of a wave often tend to maintain 
their identity as they propagate, which is represented by surfaces 
everywhere perpendicular to the dixection of wave motion. These 
surfaces are called ''surfaces of constant phase'’ or phase surfaces. 


The propagation of gravity water waves can be represented by a form 


é +> 
Cik,t) = aGt, tet Ot) (3.1) 


where a(x, t) is an amplitude function aud the sinusoidal term provides 
for the motion of the wave, where the surfaces cox , t) = constant are 
the surfaces of constant phase [Morse and Feshbach (1953), Phillips 
(1966}]. 


This physical interpretation of the phase surface function © yields 
the definition of the wave-number vector field K and the scalar 


wave-frequency field % in terms of the phase function: 


; KR = vo (3. 2) 
* , 3 and 
j = 22 ff 
L xt (3. 3) 
‘3 a In particular, the classic solution for the surface oscillation of a 
: ’ progressive water wave moving inthe +x direction [Lamb (1945), 
oe Stoker (1957), Wiegel (1965)] is given by 
Plas 
me ‘ € 
. C, (*> t) = asin an(# -+) (3. 4) 
oe q where 
a is the wave amplitude 
L is the wave length 
4 and T is the wave period. 
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Thus application of (3.2) and (3. 3) to (3.4) where p= Zn( : - & ) 


yields the wave-number inthe +x direction as 


k 


o 


tl 
N 
rls 
? 
= 


and the wave-frequency 


Ww 
o 


ill 
als 


(3. 6) 


Note that in Equations (3.5) and (3.6) a subscript o has been utilized. 
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In all following analyses this subscript refers to conditions of no wave- 
current interaction and not to deep-water conditions as is often cenoted 
in the literature. For deep-water conditions the subscript d will be 
utilized. 


Since the curl (grad @) =0, then Equation (3. 2) becomes 


vx k = 0 (3. 7) : 
and consequently the wave-number vector field is irrotational, More- 
over if p (x, t) is a contimious function then the order of differertiation 
yields identical results and consequently 


3 3 
2 (ve) = (2%) (3.8) 
Thus substituting from Equations (3. 2) and (3. 3) yiel¢s 


> 
ak =m 
3t +yu=0 (3. 9) 


Equation (3. 9) is a kinematical relationship which describes the con- 
servation of wave number. Consider a single wave train being viewed 
by an "Fulerian" observer at a stationary point. The time rate of 
change of waves viewed by the observer must be balanced by the con- 


vergence or divergence of the wave frequency w , which describes 


the flux of the number of waves. 


Consider now the case of surface waves interacting with a mean current 

> 

U. Kinematigal requirements yield that the wave frequency is given by 
wa 


s - |} => 
w=wtk-evU (3. 10) 
where the first term on the RHS is the wave: number with respect to the 


current system where 


w = w(k,x) (3. 12) 


In the following analysis concerning surface gravity waves it is 

> 
assumed that the depth of water d and meancurrent U vary slowly 
so that the classical solutions for no wave-current interaction are 


valid during interaction such that 


w? = gk tanh(kd), (3. 12) 
where g is the gravitational constant, 
the phase velocity ¢c in the local wave direction is 

c* = & tanh(kd) | (3. 13) 
and the group velocity is 


ee 2kd 
co ea a te, (1 + Sinh(2kd) 


) (3. 14) 
Figure (3,1) schematically describes the basic wave-current inter- 
action terminology. Furthermore all following analyses will assume 
that averaging over the water depth or vertical integration has taken 
place. From the condition of the irrotationality of the wave number 
vector & in horizontal space coordinates x and y due to vertical 


integration, Equation (3.7) becomes, in cartesian coordinates, 


“~ ok, ok 
Wk a aS ea = 0 (3. 15) 
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Figure 3.1%. Schematic View of Nearshore beach Terminology 
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where 


ky = keos 6 (3. 16) 
and 


ky = ksin 0 (3.17) 


Furthermore, assuming steady flow conditions exist, then 3.5 8) 
and Equation (3.9) becomes 


vw = 9, w+ K+ U) = 0 (3. 18) 


and for an arbitrary mean current system, the gradient of a scalar 
field can only be identically zero if 


> 
wt+k*U = constant (3.19) 


> 
If U= 0 then Equation (3.19) becomes identically the invariant wave 


frequency We and thus Equation (3. 19) becomes in cartesian form 


1 
[gk tanh(kd)]* + U(x, y)k cos8+ V(x, y)k sin§ = Wo (3. 20) 


where wy = 2WwT. and after substitution of Equation (3. 12). 


Expanding Equation (3.15) yields 


OU i gee meen ll ole oe 
costs, t sintay e COSR Ee Oy sin§ 


y 


k 


x (3. 21) 


RI 


where the wave number k is defined by the transcendental relation- 
ship (3.20). Notice that if a local coordinate system s and n as 


shown in Figure (3.1) are utilized, the form of Equation (3. 21) becomes 


De _ 1 Dk 
Ds iz k Dn (3. 22) 
with ox = cos6 (3. 23) 


{ 


and eI = sing (3. 24) 
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aod : where the operators of s and n are 


Ay 
1 


2 2 Os 
cos0<— + sindsy (3. 25) 


and ind = + cos8 (3. 26) 


Flo 


Eouations (3.22), (3.23) and (3.24) are very similar to the kinematical 
relationships obtained by Munk and Arthur (1951) starting from 
Fermat's principle of minimum travel time for a water wave ray or 
orthogonal except that Equation (3.22) is replaced instead by 


De _ _1 De 
Ds “cD (3. 27) 


where c is the phase speed of the water wave as given by Equation 
(3, 13). In fact, if the mean current is identically zero 

U = V =0, then it can be shown that Equation (3, 22) reduces exactly 
to (3.27, aad thus (3.22), (3.23) and (3. 24) are the general relation- 
ships governing the ray path with wave-current interaction. 


While the form of Equations (3. 22) to (3.24) appear deceptively simple s 
such that a standard numerical computational technique such as a . 
Runge-Kutta or similar method could be utilized, an expansion of the 

RUS of Equation (3.22) yields a problem. Differentiating Equation 

(3, 20) yields 


"] oe {* 32( Usine . V oss) z x( cose + sind =) i 


. 4 Pa0 
‘ - Bk sech (kd), od + {cons Vsing (3. 28) 
2[gk tanh(kd)]2 °* 


, Blkd sech"(ied) + tanh(kd)] 


2[gk tanh(kd)]* 
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and 


dk | 28 ayes aU, sine 2¥ 
ay “fe ay (Using - Vcos§) - (coro 2 + sin® ay) 


Co? 
a -&k_sech (id) od + {vcoses Veing (3. 29) 
2[gktanh(kd)]2 °Y 


, glkd sech (led) + tanh(ied 


2{gk tanh (kd) J2 


and notice that both ok and ok each have a term 96 and 98 ; 
ax oy 38x oy 


respectively. Thus Equation (3.22) does not explicitly yield an 
expression for the ray angle 6 in terms of only changes along the 
ray path s . Hence the valuable technique of integrating along 
characteristic lines is no longer valid if Equation (3.21) is to be 
fully solved. 


3.2.2 Wave Dynamics 


As indicated in the introduction, Section 3.1, the objective of this 
current research effort is to determine the effects of wave-current 
interaction on the nearshore circulation characteristics. Thus of 
prime interest with respect to wave dynamics is the change in 
wave height characteristics as the wave interacts with the nearshore 
current distribution. The conservation for mass, momentum and 
energy per unit area due to the interaction of wave motion ona 
variable current have been given by Longuet-Higgins and Stewart 
(1960, 1961) and Whitham (1962), In this section the important 
relationship arises from the energy balance of the fluctuating motion 
of a wave train in which energy dissipation is negligible. 


Vertically integrating the energy balance due to the fluctuating wave 
train superimposed on a variable current system and averaging over 


time during a wave period yields the energy relationship 
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E = z og? is the energy density per unit area (3.31) 


Oo; : is the "radiation stress'' for surface waves 
defined by Longuet-Higgins and Stewart 
(1960, 1961) 


and given by 


Cee = Elf2n-})cos*+ (n-4) sin“6] (3. 32) 
Oy, = El(2n-}) sin’@ + (n-2)cos“6] (3. 33) 
| a, ee = nsin(20) (3. 34 
xy yx 2 
where 
c 
= (8) 22 —2kd__ 
n= € ) a¢ + TanaRa) (3. 35) 


Since the region of primary concern is the nearshore coastal zone 
especially between the breaker zone and beachline, the tendency to 


consider kd << 1 as was assumed by Noda (1972,1973) is very strong 


and outwardly very reasonable. But a more careful analysis of the 


physical processes involved in the breaker zone deems this unwise. 
In particular consider the degenerate case of surface waves propa- 
gating inthe +x direction on a variable current U(x) in infinitely 

deep water. In this case since 6= 0 everywhere the kinematic 
relationship Equation (3.21) is identically satisfied and Equation 


(3. 20) yields a quadratic equation with solution 
ce 1 
- _9) 4U \2 
c= fie (1+) | (3. 36) 
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where the positive sign in the square root term is taken so that 


1 
c=) = cy -(£) (3. 37) 
° 


when U=0. Notice the interesting effect that no solution to Equation 


(3. 36) can exist if U/c,<-%. At the critical velocity of U = - 2 


the square root term becomes zero and Equation (3. 36) yields 


c= (3. 38) 


| ay 


y 
c 


and = =$ (3. 39) 


Since the local group velocity of the deep-water wave system is 


c. = te (3. 40) 


then Equation (3. 39)physically implies that the wave system can no longer 
propagate when the mean current exactly opposes the energy propa- 
gating speed of the wave system. 


For this special case the energy relationship Equation (3. 30) becomes 


RMR a a ar 
E dx (U+ c,) dx 2(0 + ¢,) dx a 
and the solution to (3. 41) is 
H c 
i (3. 42) 


o §=6 [e(e + 2U) J? 
These results were given by Longuet-Higgins and Stewart (1961). 
Extension of these concepts to the nearshore coastal zone implies 
that the complicated vector direction of the mean current coupled 


with the wave direction could yield the equivalent situation where 
the mean current directly opposes the local energy propagation 
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or group velocity of the wave. In the limit aa this situation is approached 
the local wave length will approach zero with respect to a stationary 
observer. Thus while the local water depth d may be small, the local 
wave number k = 2n/L_ may become very large such that the so called 
"shallow water" approximation may not be valid. Hence in all subse- 
quent theoretical formulations with wave-current interaction, no 
approximations are made for the magnitude of the term kd. 


Expanding Equation (3. 30) in cartesian coordinates yields 


L aE ay 
(U + c,cos8he ax + (Vt c,sin®) ay 


gO AUt c,con8) Pome: ee es c, sin®) 


ox oy 


xx Ox yx ay xy ox oyy oy 
where 
- = 1 2 1,_. 2 
O,, = (@n-2)cos 6 + (n - 2)sin'®@ (3. 44) 
= a by: 2 1 2 
Tyy = (2n-3)sin 6 + (n- 3)cos 6 (3. 45) 
- _@z _-R.; 
and Ty = " =5 sin(29) (3. 46) 


Since E is defined by Equation (3.31), substituting for E in Equa- 
tion (3, 43) provides directly a relationship for the wave height 


2. aH a) 2 8H , 3 (U+e cose) ; 

U+ c,cos8) Fs + (V+, sin@) H ay Se g i 

_ j 

+e Fee) 4 5 = 0 (3. 47) 

where : 
4 

Go = (z.3% + To5 -~ + Tay — + On ay] (3. 48) i 

: 

3 
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Finally expanding Equation (3. 47) fully yields 


2 3H { 2 3H au av 
(Urs cos8) a oa + (V+ c,sin6) oe + —— + ay 
ac ac 
- , 38 £ Pele) P g 
- sind Ss + cos8 ox + c,cosd sy + sin8 Sy 


The group velocity and wave celerity functions are given by 


i 


2 
c = E tanh ica] (3. 50) 
_¢ 2kd 
“g ~ 5} ecm] (3.51) 
oc, a k oe tad oe +] sinh (2kd) - 2kd cosh (2ke)] 
ox sinh “(2kd) 
) 2kd ac 
: [1 + sinh (2ka) 4 3x (3. 52) 
oc, e[ ae +d ak] -[ sinn (2kd) - 2kd cosh (2kd) 
= sinh “(2kd) 
ee |G eee eee 
2 sinh (2kd) } dy (3, 53) 
where 


ge. Lk 2 od, 4 ok ok 
+ [ksech (kd) k 22 + a SE) ~ tanh (kd) =| (3. 54) 


3¢ = —E— |}. 2 Oe 3d ok mi ak] 
ans [a sech (kd) (kc 5° + d 5) ~ tanh (ica) SE (3. 55) 
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and where oe and ok are given by Equations (3.28) and (3. 29) res-~ 


ox oy 
pectively and k defined by the solution to Equation (3. 20). 


To understand some of the physical processes interacting within the 
energy equation (3.49), it is useful to transform (3.49) in terms of the 
local coordinate system s along the wave ray and n along the wave 
front. Utilizing operators defined by Equations (3.25) and (3. 26), 
Equation (3. 49) becomes 


1 
H Ds Dn ec, Ds ec, y oo ox ay 
a 
+ 2c, =0 (3. 56) 


If the mean current velocity is now set equal to zero, U=V=0, 


then Equation (3.56) becomes 


De 
J1DH,1D0, 1 _‘g. 
H Ds *2 Da* 2c) De 7? (3. 57) 


or in terms of the energy density 
2 DE, DO, 1 Mg 24 (3.58) 
& 
Notice that both Equations (3.57) and (3. 58) describe the changes in 
the wave height or energy as being related to the curvature of the 
wave front Be and to the logarithmic change in group velocity along 
the ray path, = and - In point-of-fact ne describes ray refrac- 


c Ds 
4 
1 De 
tion and =. 5° wave shoaling. 
g 
The form of Equation (3. 57) suggests a form of separation of variables 
where = } 
H = H,(6)H ,,{¢,) (3. 59) 
and substitution of (3. 59) into (3.57) yields 
DH l De DH 
tt —_£2 1 wk, 1 DO. 
HDs ‘2 De +H. Ds * 2Da~° 13580) 


ted onalb. dl idian ttle. da aul 


Pa the ted 


Bobeirtiedteteaciteo cit. qcegdbltMins penttitlea’ series eo eer ee 


et Deel a aS tg 


Since Aon is only a function of c 


, and H,. only a function of 9 , 
Equation (3,60) implies that 


1 DH,» l De 
en ——- a8. = 1 
H Ds + 2e Ds c (3. 6!) 
sh 
DH 
l r 1Dd6 _ 
and H_ Ds *2Da~ -° 


(3. 62) 
where C is a constant. Equation (3.61) can easily be integrated and 
yields a solution 


— 


c 
Hon = (3. 63) 
c 
j g 
\ ‘ 
, 4 where Cc, is a new constant. The solution to Equation (3.62) is much 
\ more complicated. By considering the ray separation diagram shown 
in Figure 3.2, Munk and Arthur (1951) have shown that 
\ 
\ og 1 Db _ DO 
\ b Ds = Da (3. 64) 
\ 
\ and defining £ = b/b a? the equation for ray separation (3.65) becomes 
7 1 DB _ Da 
4 ) Ds Dna (3. 65) 
Now substituting Equation (3.66) into (3. 62) yields 
DH. 
at By DB 
H, Ds * 2p Ds ~ “© (3, 66) 
and integrating directly yields 
4 Cc 
4 H, = mo (3. 67) 
8 


where c, is a constant. 


Thus finally substituting back into Equation (3. 59) produces 
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Figure 3.2: Ray and Wave Front Terminology 


H = (3. 68) 


and in deep water H=H, ’ Cyr c and B?1 yields 


bq 
C=H ¢ (3. 69) 
oO bq 
and finally, S 1 : 
72 ( =a)" (3. 70) 
ONeg sh YB 


where Equations (3. 50) and (3.51) give 


1 1 


“Ba\ ? 1 7 
=*) = (3. 71) 
( “g ) ae h(kd) f + shor st 


Equation (3. 70) is the well known classic solution to waves undergoing 
transformation due to both shoaling and refraction, 


The solution for H from Equation (3, 70) is not fully complete since 
B is yet an unknown. Munk and Arthur (1951) have derived a dif- 


ferential equation for 8 , called the wave intensity 


2 " 
2B + p(s) BE + qieip= 0 (3. 72) 
Ds 
where 
l P 
p(s) = -cose[ = 8¢] - sine| 5 $e (3. 73) 


2 2 2 
and q(s) = vin’ of 7 25] - 2 sindcose| ¢ 5 + cos“a[z al 


and solutions for 8 are shown in Noda (1972,1973), In particular 


it can be shown that Equation (3. 72) degererates to the Snell's Law 
solution when d=d(x) only 
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i 
_jcos® 2 
8 les a (3. 75) 


Returning back to the general equation (3.56), it is painfully evident 
that with wave-current interaction, the simple concept that the local 
wave height can be represented by a product of wave shoaling and 
wave refraction factors as described by Equation (3.59) is no longer 
valid. The combined dependency of refraction on shoaling and vice 
versa necessitates the solution of Equation (3. 49) directly. 


As the wave propagates from relatively deep water into shallow water 
or into an area where mean current conditions exist, Equation (3. 49) 
will govern the local wave height until an instability occurs. This 
instability is usually wave breaking due to the effects of shoaling, 
refraction and wave-current interaction, In order to determine when 
breaking occurs it is assumed that spatial variation in U and V 

are sufficiently gradual so that an empirical breaking criteria is 
imposed, developed from the non wave-current interaction observation. 
The theoretical limiting wave steepness condition from Miche (1944) is 


. TT 
" = 0. 142 tanh (=) (3. 76) 


where the subscript b indicates breaking conditions and the breaking 


wave length L, is given by 


= 20 
Ly *t, (3. 77) 


and k, is derived for the transcendental Equation (3.20). An examin- 
ation of experimental data of waves breaking over a horizontal bottom 
by Le Mehaute and Koh (1967) indicates a better limiting steepness 
criterion is 


A 


zo = 0.12 tanh ( =) (3. 78) 


By 
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Since the wave number k during wave-current interaction is obtained 


directly from Equation (3. 20), then Equation (3. 78) can be transformed 


“ H, = a tanh (72) (3. 79) 


; During computation if the solution for the local wave hcight H from 
_ Equatior. (3. 49) is less than Hy; then the wave height is H. If 
computation indicates that 


H 2 H, (3. 80) 


then the local wave is considered to have broken and the empirical 
relationship Equation (3. 79) is imposed where H = H- Thus the 
effects of the mean current become critically important through the 
wave number k. In other words if the local mean current is in the 
same direction as the local wave direction then the local wave number 
k becomes smaller which requires a larger wave height for breaking 
to occur. On the other hand if the local mean current opposes the 
local wave direction then the local wave number k increases and 

the limiting local breaker height decreases, This phenomenon is 
easily seen at the entrances of river and estuaries when an outflowing 
current meets an incoming gravity wave system. The local limiting 
breaking wave height decreases so that even very small waves seem 
to "white cap" and break, 


Application of the empirical breaking criterion Equation (3. 79) is 
indeed crude. Recent studies of breaking waves by Divoky, Le Mehaute 
and Lin (1970) indicate that wave breaking is dependent on a ''charac- 
teristic" bottom slope and the research effort of Galvin (1969) centers 
on different types of breaker characteristics as a function of bottom 
slope. Thus the breaking criterion expressed by Equation (3. 79} 
hopefully indicates the major breaking process, sufficient to deter- 
mine the merits of the concept of nearshore wave-induced circulation 


including wave-current interaction, 


Ae eon ce ee 


3.2.3 Nearshore Wave-Current Circulation Forrnuiation 


In the following theoretical formulation, the equations which describe 
the circulation pattern within the nearshore zone are derived. The 
formulation objective is to initially solve for the wave height and 
direction fields and subsequently for the wave-induced nearshore 
circulation pattern with no wave-current interaction [Noda (1972, 
1973)]. The next step is to use this circulation pattern or more 
specifically the circulation velocity fields U(x,y), V(x,y) as the 
mean current to be input into a new calculation of the wave height 
and direction fields and again obtain the solution for the new wave- 
induced circuirtion pattern with wave-current interaction. In theory 
this iterative technique, which assumes a series of quasi-steady 
states, should hopefully in the limit, converge such that the circula- 
tion velocities are exactly equal to the previously imposed mean 


current. 


The coordinate system is described in Figure 3.1. Vertically inte- 
grating the momentum and continuity equations and neglecting the 
nonlinear and time dependent terms yields 


1 oo, eT. 
Sax“ -ptd| ox * =| - F, (3. 81) 


an 1 30 oT se 
& oy = > pint a) “55° + = - (3. 82) 


and 2 [ulnt+d)] + = [vint+d)] = 0 (3. 83) 


where 
7 is the mean water surface 
e is the fluid density 


Fis y are bottom friction terms 


and u,v are the new wave-induced circulation velocities, 
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The analysis of the bottom friction term is similar to that provided by 
Longuet-Higgins (1970) and a detailed formulation of this concept is 
given by Noda (1972) where 


FL Pre (3. 84) 
and 
Fy REX vaioe ce (3. 85) 
where 
G is a constant bottom friction coefficient, 
usually ¢ = 0.01, 
and T is the local wave period. 


At suming that 


ntd d 


and cross-differentiating Equations (3, 81) and (3. 82) and introducing 


a stream function 4 defined by 


ob =. 
and 
ot = t+vd (3. 88) 


automatically satisfies the continuity equation (3. 83) and yields the 


nearshore circulation equation 


ai ,%y , LaF au, 1 aFay 
ax” ay” F ay oy F ox 3x — 
taf do, 30 5 T1232 aa 
F 2-13( 8x sf oy ) ° 2 af ay . ox )] 
where (3. 89) 
Fo = ch (3. 90) 
aT sinh kd 
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2 
of” = LE (n- 4) cos%e + (n-3)sin?9} (3.91) 
Z 
¥ H 3 
oy, = B-[(2n- 2 sin’@ + (n-$)cos"9] (3. 92) 
S ASe = are bd 2 
icy = Tyne 7 Te Bein(28) (3. 93) 


Typically since the nearshore coastal zone is of primary interest, the 
urge to utilize the so called ''shallow water" approximation where the 
argument kd is considered small seems appropriate. Note the 
important concept distinction that vertically integrating the momentum 
and continuity equations is not identically synonymous to the shallow 
water approximation. As described in the previous section, while the 
water depth d may become small, the wave number k may become 
very large if the mean current opposes the wave ray direction and its 
magnitude approaches the energy propagating speed of the local wave; 
implying that the product kd may become very large. Thus the form 
of the circulation equation (3. 84) should definitely consider this inter- 
active concept. In particular the local wave period T in Equation 

(3. 90) is not invariant to a stationary ''Eulerian" observer and its 


variation must be considered. 
Since the local wave period T is defined by 


= 20 
T=2, (3. 94) 


then substituting Equation (3, 94) into (3. 90) after noting that w is 
defined by Equation (3.12) yields 


p - SV Hal (3.95) 


ac sinh (2kd) 


and consequently its x and y derivatives are 


ti) 


SNe te A AC ta pS FEE Ra EES RA LR BE EDS REE RS ES TESTS: AR NTO RET A G/ERESE 27 ERS AE 


i 
Zr si ok od dk 
QF _ tV2e H Lae veal {_ oie 2k\ k aa + ds= | 
T 2 k 
2d 


ox sinh“(2kd) 
1 
ZPqch _ op od 
- k dsx 2H ox (3. 96) 
sinh(2kd) 2 . 
£5 ak ad, 4 3k 
aF _ tV2e H_[ sinh(2ka) Be sinh(2kd) ae 2ki k 5 +d 5 cosh(2kd) 
oy T U2a2 k sinh“ (2kd) 
zfach ~ 23d 
,7/—_ k_ _oy oy (3. 97) 
= sinh(2kd) 3 ' 
and thus 
( ad ak } 
ter ade OS oe Ses) pdf at gta: ~ pa 5a 
F dx ~ 2k| dx tanh(2kd) Hd\¢ 3x 3x ; 


od ok 
dar _ ipa May * Cay) ) 1 (48H on 24) 3,99) 
F dy ~ 2k] dy tanh(2kd) Hd\° ay oy ‘ 


where k is the solution to Equation (3, 20). 


Carrying through the derivatives of the RHS of the nearshore circula- 


tion equation (3. 89) yields 


2 42 32 32," 
1 xX x 
RHS [Eq. (3. 84)] = &{ = |—* + — 3 . —_ YY. 
oa 30° an ed aan ac 
ee) xX 4 XY), OX} vy , xy 
a’ Ox oy ae ay Ox 


(3. 100) 
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and full computation of the derivatives in Equation (3. 100) yields 


ea 2 

xx _H 38 on 2 2 3H * 

Sa OB: [-»38 sin(26) + 35 (1 + cos al tT Sa ax (3. 101) 
a0 He 


n $2 sin(26) + 3a (1 + vin”e)| + wo (3. 102) 


ok 
oT 2 we 
et E [ 22 cos(29) + & n SE gin cea] + Lon T 


Ox eo) nox "yx (3. 103) 
an 2 
H a6 n 0H lon ,* 
23) gee ope 
sy S| 2 3y cos(26) H dy sinize)]+ 2 a oy T ae (3. 104) 
where 
ad ok 
dn (: k 3n * Tax «) 1. —2kd (3, 105) 
x setakah tanh(2kd) : 
on ko + 2ked 
—= l- ——yprT (3. 106) 
oy sinh(2kd tanh(2kd) 


and the second-order derivatives are 


2 * 
a0 2 2 
be eee 2 ae ee pe _ an 28 38 _ 36 dng; 
8 n sin(26) - 2n ay cos (24) Sx dy sin(2§) 


oyox ayex Ox 
98 on a: 3p 
Fay ox sin(26) + (1 + cus’? > al 
+ HOH) 88 cincae) + ito 2g) on 
4 y Ox Fs \ US 5x 
ac. 2 
2 dH 0° xx + [i *H 1 aH 3H 
+ Ff Ox oy e206 E oyax — He °%* | (3-307) 
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2% 
ate 2 2 
XY - 2)? [23 cos(26) - 2(22) sn + 4H 30 2H cos(20) 


2 2 
+ [= as +(2#) ssc B 22 % cos(2@) + Hi envo] 


ay y ay 
ar 2 2 
+ Ldn wy, Ps 133+ i) (3. 108) 
n oy xyin dy n y 
v6, H“}  aga0 ae 36 an 
tt = Fn $228 cosi2a) + Se ain(20)] + S> SP sin (26) 
, an 28 an 
oy ax a 


+ ke sin(26) + 32 (I + sin ‘a 


ao" ; 2 
2 oH My ,,¢ [1H 1 aHaH 
* Hay ox 34 2a | oxdy ~ 42 dy 2 rae 
arr 2 2 
: 
—e . 3/2 ¥ cos(28) - 2(88) since] 4H $8 2 ens (26) 
ox ox 
a¢H . faH\* H an 
+1H sue (2) sin(29)> + 3 Sx | Ee 88 8 cos(28) 72 es H sin(2@ 
x. 
xr 2 2 
its t [hee 4 (zy 110) 
ae ee ee YX" ax n wy 
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SiR gS NGA SNL OSG EGP GTS PETES LET PGE IED CLAD TEAS TOE INS | 


ao 2 x 
0c H CT] 2 aH 
= -n 5, sin(2e) + (1 + cos 29) 22 +2 9H: (3. LEL) 
oy S 2]. H by a 
and 
a6 2 
* 

ae ae E 98 sin(26) + (1 + sin?6) a2 +238 goo (3. 112) 


Notice all derivatives and parameters have been cpeci ee rs for 


2 
second-order derivatives in n , i.e. on : 2n and wn . 
3x7 By" axdy 


While algebraically feasible, these derivatives would contain 2nd 
order derivatives of k with respect to x and y and considering 
the complicated form of Equation (3. 28) and (3. 29), it was decided 
to use central differences of the first-order derivatives to compute 
the second-order derivatives of n. 


Finally the RHS of Equation (3. 89) and the variable coefficients of the 
first-order derivatives of 4 are completely specified such that 
providing sufficient and necessary boundary conditions, the stream 


function can be found by utilizing an iterative relaxation technique. 
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3.3 Numerical Analysis 


This section deals with the numerical techniques developed to deter- 
mine, first, the wave characteristics of direction and height due to 
wave-mean current interaction,and second, to utilize these results 
to determine the resulting wave-current induced circulation pattern 
in the nearshore coastal zone. Hence, essentially two separate pro- 
grams exist--one to obtain the wave local height and direction, and 
the second to solve for the stream function 4 once the wave height 


and direction fields are known. 


In a previous study [Noda (1972, 1973)], the solution to wave-induced 
nearshore circulation was found excluding the effects of wave-current 
interaction. In that study, the wave height and direction fields were 
obtained by integrating along characteristic lines by utilizing a fourth- 
order Runge-Kutta scheme. As discussed previously, with wave- 
current interaction, characteristic lines no longer exist and thus a 
completely new numerical technique was developed. This technique 
first solves the kinematics problem directly on numerical grid points 
i, j and yields the ray directional field @ at all nodal points. Then 
the energy equation (3. 49) is solved directly on the same numerical 
grid points i, j to obtain the wave height field. 


The numerical technique discussed above to determine the wave char- 
acteristics for wave-current interaction makes important use of the 
fact that prototype data indicates a longshore periodicity of both the 
wave -induced circulation pattern and the bottom bathymetry such 
that 

d(x,y) = d(x,y +m) where |m| = 1, 2, 3 (3. 113) 


This longshore periodicity leads to the key boundary condition which 
emits the development of a highly efficient numerical algorithm for 
the solution of wave characteristics including wave-current interaction 


in the nearshore zone. 
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While longshore periodicity proved to be a valuable tool for the work 
herein, nevertheless, the method of integrating along characteristic 
lines is such a generally powerful method that continuing research 
efforts should be extended to be able to utilize this technique for 


wave-current interaction, even as an approximate approach. 


3.3.1 Numerical Solution For Wave Characteristics With Wave- 


Current Interaction on a Longshore Periodic Beach 


Figure 3.3 schematically describes a longshore periodic beach with 
wavelength 4}. Equation 3.21 describes the variation of the wave 
direction field 6 asa function of x, y, U, V, and T where the 
wave number k is defined by Eq. 3.20. Rewrite Eq. 3.21 in the 
form 


38 1 Bk) y ang [284 2 Bk] - 
cos 0 [28 ae se | + gin 8 EE + kK ak | = 0 (3. 114) 


and rewrite the derivatives of K from Eqs. 3.28 and 3. 29 as: 


1 a _ 36 [U sin@ - V cos@] 1 6k 
k 6x 0x A * kk ax ieee) 
1 dk _ 06 [U sin 6 - V cos 6] i ok 
Ee * 3p = Fs (3. 116) 


1 2kd t, 
where A = Ucos@+ Vsin@ +> [a+ ginh(2kd) | [ =2-ucose-Veine | 


(3. 117) 

Ede =~ [3 coe Fy sine | -[7o- uncoss-vising ]38 
sinh(2kd 

A (3. 118) 


r3u av .. 


sinh (2kd) 


A 
68 (3. 119) 


BEACH LINE 


OFFSHORE 
“BOTTOM 
oy ™~. CONTOURS 
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@ WAVE CREST 


: Figure 3.3: Schematic Iluatration of Periodic Beach 
ro Terminology 
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A -2-5795 


F Now Eq. 3. 114 can be written: 


36 sin 6 (U sin§@ -Vcos@ 30 /[.. cos @(U sin@-Vcos8) | _ 
se [coa0+ A sos) | +B Laine - A ]- 


fou fey 
oF 


l 1 dk. 
k yee % dx sin 8 


(3. 120) 


Thus, at any grid point i, j utilize a forward difference derivative 


in x and a backward difference derivative in y and solving fer @ 


yields , 
—= —_ Aaa Q. 
abo of Take 13k Ee beats FS 
8 By; Kay cos 8557 k ox sin Q; ;+ sin®. .- RS 
Orn mite 2 7 
(U sind, - Vcos 0. iD LL cos 6. ij U sind. ig ¥ 008% ra 
(3. 121) 
Bere : sin 6; P aca (U sin 8, Tei —ne ze hi] 
ij Ay 
(3. se 
and i and + = are evaluated at i,j. 
The RHS of Eq. 3. 121 contains terms cos, and sin Q j and a 


variety of approximations for these terms can be utilized to update 


8; iG The simplest would be to use the previous value of 96 A 
more sophisticated approach is to approximate these sinusoidal 


function at i,j in terms of the four surrounding grid point and 


going to 2nd order using Taylor series expansions yields 
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P _1f.. ‘ . : 
sinQ; 5 = Zz aad Fae ainQi sing 4,t eine, | 


1 
+e essa g 4-1,9) aa oa cc 


+ (8; 544° 94,51) (Cos 8. |- cos aad) (3. 123) 


me 7 
cos 6. aay [cos 8543 3+ cos ae cos 0; jtat cos 8s j-1d 


+ 


l A : 
#1 5 ie ee (sin O44 jo sin 85 -1,)) 


+5 5417 41,7) (29M jay sine, y) | (3. 124) 


e- 


The need for the sophistication of Eqs. 3. 123 and 3. 124 is question- 


able, although utilizing this scheme to iterate for 0; ; yields amazing 


, 


results. For instance, if d= d(x) = mx (a plane beach), then starting 
with a boundary condition cifshore and setting the initial 6 field equal 
to m, a solution using Eq. 3. 121 will conveize to within 1% of the 
Snell's Law solution after only two iterations through the full field ! 
Note that it is not necessary to relax over the whole field, but the 
solution could have been obtained by iterating row by row or along 

i = constant , Figure 3.4, working inward toward the beach. In this 


c.se, the approximation given in Eqs. 3. 123 and 3.124 for sin cP j 


and cos 95 j should not contain values of § on rowi-l. Neverthe- 
less, ‘he extremely rapid convagence of this technique is suffici- 
ent (o justify the numerical algorithy. The row by row technique 
wi.l be utilized to find the wave height field, shown later in this 


section. 
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Figure 3.4° Local Grid Description 


Before further exploring the numerical techniques, the important 
longshore boundary conditions should be examined. Figure 3.5 is 
a grid description of the entire field of calculation. Notice that 
while the beachis periodicin \ orfrom j=2 to j=rtl,a 
column has been added on either side. Thus, computations for 
either 9 , H or Wy , are only performed inside of the boundary 
lines from i= 2 toi=m-1l and j=2 to j=rt+1. To utilize 
the important longshore periodic boundary condition, for any com- 
puted variable £ along each row i the following conditions are 
imposed as the computation proceeds toward the shore line in 


decreasing values of i: 


f +17 4,2 (3. 125) 
i eae “ys (3. 126) 
t l = fF m (3. 127) 

Imposing conditions 3,125 to 3.127 upon 6, Hor y physically 


requires that the functions and moreover its derivatives are continu- 
ous and periodic in \. Hence, iteration continues for 8. j until 


every updated vaiue of 8; j is less than a specified percent of 0; j 


itself. A typical run will converge in between 3 to 7 iterations with 
a maximum relative error for each a. . of 0. 1% of itself. 
At this point, it should be noted that for each updated value of 0; j 


a transcendental relationship for k must be solved defined by 
Eq. 3.20. To efficiently solve for k the Newton-Raphson method 


wag utilized where 


2 
= T + 
e(k) = gk tanh (kd, ) - [ - Ukcos@, ; - Viesin a5] (3. 128) 


1 his eat NoPE TT ge gm eS Cag 
ail dln hic eh i 


SHORELINE 
Full Grid Description 


x 
Figure 3. 5° 
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and 
de = 6! (k) = ¢ [kay sech” (kd) + tanh (kd; | 

+2 [u cos8, + V sind, Ji -Ukcos 8, .-Vksin8, Bi (3, 129) 
| rj i,j 

Hence, (x ) 
Mie s koa rs : ic } (3, 130) 

old 

and this iteration was performed until 

Le - Koa! = 9.0001 [« ew! (3. 131) 


Computation starts far offshore where the periodic beach d(x, y) is 
defined to be a plane beach d = d(x) = ax starting at i=m and 
offshore from this row. On the row i-=m, the local wave angle 

is specified in reference to a chosen deep water wave angle, using 
Sneli's Law. From this point, computation immediately begins with 
Eq. 3.121 and the output is the direction rar for a given relative 


accuracy. 


oad The next series of calculations solves for the wave height field. 
: Rewrite Eq. 3. 49 to yield 
oH _ oH » 4 08 30 
+ =>: — 
a a: + {Vv Ry 7) 


(3. 132) 


Define 


28 cos 28 au, 2) 


ax “g dy ‘tx ay 
oc ac 
- cos = - singe - a} (3. 133) 
where (38) and (22) are obtained by central differences, 
i,j i,j 


and taking a forward difference derivative in x and a backward 


difference in y and solving for Hi ij yields: 


(V+c_ sin 6) | - (U +cg cos Q) Aiea; 


a Ay Ax 
[2c vig aie licen Oo (3. 134 
i,j (Vte gin 8) ; (U + ¢_ cos 6) an 
Ay Ax 2 


In the computation for the wave height Hy i’ a row by row relaxa- 


tion technique is utilized starting on row i =m-1, and proceeding 
inward to row i=2. On each row, the boundary conditions 3. 125 


and 3.127 are utilized and the solution is reached when 


IH ow ~ Hog! £ (0-000) em (3. 135) 


The convergence of this scheme is amazingly rapid for even the most 
complicated bottom bathymetry and mean velocity distribution. Usu- 
ally, only 2 to 3 row iterations are required to meet the criterion 
defined by Eq. 3. 135. 


Similar to the §@ calculations, the wave height calculations start 
at row i=m-l where values of the wave height on rowi=m are 
used, derived from the Snell's Law relationship. During each update 


calculation of Hy je the breaking height i+ also calculated, and if 


76 


H, j exceeds this value, then the wave breaking height is automati- 


s 


cally imposed and a flag is set to denote if the breaking height condi - 


tion is being utilized at the point i,j. Computer outputs of this 
breaking index will be shown later, 


Now that the H and 6 field has been specified, the numerical 
algorithm proceeds to the computation of the stream function field h. 
The technique used to find the stream function is very similar to the 
technique used by [Noda (1972, 1973)]. Equation 3. 89 is rewritten 


2 2 
3 3 loFr\sa 19a fe) 
oo + So + (GE) M+ (GM. w (2.430) 
: x ay 
where x a2q* a2, = 
7 ae on 2. _XY_xy. 
“*F £1d [ Syae a a " @xey xe ] 
ad d 
ay | Eee op! May le oe | Oop Ras at 
ae SX oy a2 3 


(3, 137) 


Utilizing a central difference form for the ¥ derivatives, and 


gathering terms of J. yields: 


(ay)? ZAy 
(; 2F) 
2 Fa 2 ae 
+ (ad 4) fd vi, 
th a (3. 138) 
Lay? ZAy ] 
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Note that before computation of § begins, the n field and its deriva- 
tives are first computed. Thus, the second-order derivatives for n 
can then be calculated from the first-order derivatives using central 


differences. 


The longshore boundary conditions for | are given by Eqs. 3. 125 and 
3.127, and at the beachline | is defined to be } =O. The final 
boundary condition is to move the final offshore grid row sufficiently 
far from the nearshore zone so that its influence on the nearshore 
circulation pattern is small. At this finaloffshore boundary i = mm, 
the stream function is again defined to be | = 0. With these boundary 
condition the iteration for can begin using Eq. 3. 138, and the crier- 


ion of convergence is assumed when 


|< 0.001] ¥ “3, 138) 


ly 


new ” Yoid nel 
In comparison to the numerical techniques for 6 and H, the con- 
vergence of the field requires more extensive iterating, usually 
300 to 400 complete field relaxations before the solution converges 
as defined by Eq. 3. 138. Note that condition 3. 138 is not a com- 
pletely sufficient condition to assure convergence since the solution 
may be asymptotically convergent or perhaps even divergent. There- 
fore, value of the relative error are varied to insure that Eq. 3. 138 


does indeed provide an acceptable convergence criterion. 


In the zone between i=m andi=mm, a plane beach profiie is 
assumed so that calculation for 9 and H are explicitly determined 
from the Snell's Law relationships once deep-water wave character- 


istics are defined. 
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3.3.2 Verification of the Wave-Current Interaction Algorithms 


To test the numerical algorithm for the determination of the wave 
characteristics, the degenerate cese of a wave propagating in the -x 
direction in deep water over a variable current system U(x) was 
utilized. The solution to this case was first given by Longuet-Higgins 
and Stewart (1961). For the case of a wave propagating in the -x 
direction, 6 = 180° and the analytic solution for the wave celerity 


is: 
“9 4U 
e=2[1+a-42)] (3. 139) 
Oo 
where T 2 
CA SOs iE 
co = qe = (L£) (3. 140) 
oO Le) 
- 2 
but since k, = ~o (3. 141) 
4 
g . Bo 
c= 7 a (3. 142) 


The wave height is given by: 


Hic 
(opate) 


H= - (3. 143) 
[ele - 2u)}!/2 


A simple form for the mean current velocity was chosen to be 


U (x) = 0.2x+b (3. 144) 
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ie 7 


t 
| 
| 
{ 


Essentially, two verification tests were performed. First, the 
energy equation 3. 49 was degenerated into an ordinary differential 


equation given by: 


dc 
[- $0 +) 
dH  H 2 dx dx 
dx ~ 27 =e, on 


and solved by a 4th order Runge-Kutta technique. Second, the finite 
difference Eq. 3. 134 was utilized directly to obtain a solution for H. 
The final comparison was to check each solution directly with the 
analytic solution Eq. 3.143. Table 3. 1 shows the analytic solutions 
for the wave celerity c and wave height H from Eqs. 3. 139 and 

3. 143, respectively, for Ty = 4 seconds. 


An arbitrary offshore point was chosen so that U(x) =0 and H, =1.0 
and as the solution proceeded shoreward both the Runge-Kutta techni- 
que and the finite difference scheme did not yield results comparable 
with the analytic solution. Table 3.2 describes the numerical solu- 
tions. While the difference between the Runge-Kutta solution with 
step size dx = -0.lm and the finite difference scheme may possibly 
be acceptable, comparison to Table 3.1, the analytic solution was 


unacceptable. The problem is due to the unrealistic starting condi- 


tions imposed at U(x) = 0 where ote is assumed to be zero. But 
there is an obvious discontinuity in the derivative of wu at the point 
where U(x) =0 if ae is set equalto zero. Thus, a second series 


of calculations were performed at a start point where U(x) = -1.0 
m/sec. The Runge-Kutta program then computed its own starting 
derivative and the starting wave height was obtained from the analy- 
tic solution H tart = 0. 77475 meters. Table 3. 3 describes the 
results. Comparison of the Runge-Kutta sclution with the analytic 
results of Table 3. 1 show an exact correlation and the results of the 
finite difference solution compare very closely with the Runge -Kutta 
solutions. Appendix A contains the Runge-Kutta computer program 


for the above case. 
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TASLE 3.1 


Analytic Solution For Wave-Current Interaction: Check Case 
Ty = 4 seconds, g = 9.80621, k,=0.25161618, oc, = 6. 2428272 


H/H c (m/sec) 


: U (m/sec) 


Eq. (3. 143) Eq. (3.139) 


1. 615187 
3. 761628 
6. 417806 
13. 424308 
1.0 

0. 774750 
0. 648222 
0. 564575 
0. 504132 
0. 457911 
0. 421147 


4, 992348 
3. 737028 
3. 379949 
3. 187839 
6. 242827 
7. 119669 
7, 836164 
8. 457301 
9.013317 
9. 521207 
9.991652 
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TABLE 3.2 


Numerical Results for the Wave Height H 
For the Test Case 


dHo 


Mean Runga-Kutta Solution = a 0.0 


Current 
U (m/sec) 


Finite Difference 
Solution 


Grid Size 
dx = 0.01m 


Step Size 
dx = -0.Im 


Step-Size 
dx = -1.0m 


1. 1.0 
0.99471 0. 99365 
0. 98848 0. 98742 
0. 98232 0.98129 

0. 78302 0. 

0. 63514 0. 

0. 57060 0. 

0. 50952 0. 


TABLE 3. 3 


For the Test Case 


dHo 
dx 


Mean 
Current 
U (rn/sec) 


Runga-Kutta Solution = 0.0 


Step-Size Step-Size 
dx = -0. 10m dx =0.01m 


Numerical Results for the Wave Height H 


Finite Difference 


Grid Size 
dx = 0.005m 


Solution 


0. 77475 
0. 77156 
0. 76840 
0. 76527 
0. 64822 
0. 56457 
0. 50413 
0.45791 


0. 77475 
0. 77156 
0. 76840 
0. 76528 


0. 77475 
0. 77156 


3.3.3 Numerical Results for Wave-Current Interaction 


The periodic bottom bathymetry used to test the influence of wave- 
current interaction is given by: 


exe 


d(x, y) = mx E + ae aint? (y - x tan a) | (3. 146) 


where the first-order derivatives are given by: 


x /*) 
od _,, . Ommax ynae : sin? -x tan a) cos =( -x tang) 
ox x rn ay se 
1/3 
-X 
(=— 1/3 
b ) 10 
+ mae sin y (y-x tan a) — | (3. 147) 
and <) 
b 
ge = Onmax e sin? © (y-x tan a) cos ¢ (y-x tan a) 


(3. 148) 


where the constants are given by: 


m = 0.025, a= 20 meters, % = 80 meters 
1/3 ta 
a = 30°, b= {20} wintepe !- 


The computation starts by first assuming no wave-current interaction 
exists such thatinitially U=V-=0. The wave height H and 4 
fields are obtained and the stream function y solved for. The algor- 
ithm then computes the circulation velocities defined by Eqs. 3.87 and 
3. 88, using central differences. These computed circulation velccities 
are now the mean current system which must now interact with the ori- 


ginal incoming wave system. 
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Attempts to directly impose this derived mean current system in au 
interaction process with the incoming waves leads to failures of the 
technique because the mean current system derived for no wave- 
current interaction is too large. .fence, conditions arise where the 
local wav2 is no longer able to propagate into sorme areas and k 


passes through zero and becomes negative. 


It is important to understand the specific processee involved in the 
numerical algorithm. In particular. under prototype conditions as 

a constant wave system begins to attack a coastal area, an instantan- 
eous interplay of wave characteristics, bottom sediment movement 
and wave-generated current effects exist simultaneously until some 
type of equilibrium" condition exists where the dynamic and kine- 
matic requirements are all satisfied. The numerical model does 

not aliow changes in bottom con‘iguation. Moreover, the nearshore 
circulation system may be sensitive to large changes in the wave 
height and direction field due to the instantaneous application of the 
fuli mean current system derived from the noninteractive case. 
Thus, attempis were made to proceed much more slowly by multi- 
plying the noninteractive current field by a constant less than cue. 
This still preserved the contiruity conditions of the mean flow, but 
yet allowed interaction to take place. This technique then envisioned 
some type of step by step series ui quasi-steady circulation solutions 
until the full current system could be applied and the interactive 


results yields the input current system. 


Unfortunately, while seemingly a logical procedure, the wave-current 
interactive system is sensitive to the rate at which the mean current 

is applied. Presently, the maximum interactive current applicable 

to wave-current interaction is about 50% of the noninteractive circul- 


ation system. 
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Figure 3.6 graphically shows the stream function solution y for the 
no wrve-current interaction case where Hy = 1,0 meters, Ty = 4 sec., 
a4 = 150°, and the depthis given by Eq. 3. 146. Tables 3. 4 through 

3. 8 provide the wave direction § , wave height H, breaking index IB, 
u velocity and v velocity field for this case. Notice the existence of 
the counter eddy field defined where $< 0 , which is the degeneration 
of the normal incidence negative 4 field. Alsc, all above and follow- 
ing tables will only include the region of data where Eq. 3. 146 applies. 
Offshore from this row at x = 200meters, the Snell's Law region is 
utilized and the H and @ field are easily derived. For ail cases 
herein, the maximum distance offshore is x = 345 meters. The maxi- 


mum meandering "rip-current" velocity is about 2.7 meters/sec. 


Figure 3. 7is the solution for the stream function field * with wave- 
current interaction where only 20% of the noninteractive current velo- 
cities from the solution in Figure 3. 6 were utilized. These noninter - 
active velocities are shown in Tables 3. 7 and 3.8. Notice that the 
counter eddy field has decreased a little but the general pattern of the 
circuiation pattern remains much the same as the noninteractive case. 
Tables 3.9 through 3, 13 provide the values of 6, H, IB, u and v 
for this case. For this case, an examination of Tables 3. 12 and 3. 13 
shows that the maximum meandering rip-current velocity has been 


reduced slightiy to a maximum value of about 2.5 meters/sec. 


Figure 3.3 is the solution when for the stream function field with 
wave -current interaction when 50% of the noninteractive current 
velocity field from Tables 3. 7 and 3.8 are utilized directly. Tables 
3. 14 through 3. 18 provide the resulting data for the spacial fields 
0, H, IB, u and v. 


86 


ee a 
A eee eee Ae 


+ 


Cee ae ee 


Pn, See ee 


The resulting stream function pattern from Figure 3. 8 is indeed 
startling and unexpected. The counter eddy field has now become 
stronger, and while the outgoing meandering rip current is similar 
to the previous case, a very strong inflowing meandering current 
now exists. The maximum magnitude of this in rip-current is about 
4.1 meter/sec. Presently, it is unclear what is causing this sudden 
change in the circulation pattern, since the 20% interaction case 


produces no significant effects. 


One possible influence could be the sensitivity of the circulation 
pattern to the irmposed bottom topography. In reality, the interaction 
of the ability of the bottom to change formin accord with the current 
intensity may be as critical a factor of consideration when the full 
wave-current interaction problem is considered. Another possibility 
is that the stream function solution is sensitive to the magnitude of 


the mean current and too intense a current produces spurious results. 


Figure 3.9 shows the stream function solution ¢ when the wave 
height Ha has been reduced to H 4-0. 5 meters for no wave-current 
interaction. The circulation pattern appears reasonable with respect 
to the bottom contours with a maximum outflowing velocity of about 
1,6 meters/sec. Tables 3. 19 through 3. 23 describe the resulting 


spacial variables 6, H, IB, u, and v field respectively. 


Figure 3. 10 shows the stream function solution t with a wave-current 
interaction of 50% of the noninteraction case described in Figure 3. 9. 
For this case, 50% of the velocity field shown in Tables 3.22 and 3. 23 
were interacted with the original incoming wave system and the results 
for the variables 6, H, IB, u, and v fields are given in Tables 
3.24 to 3,28, respectively. Again as in Figure 3.8, the results show 
extreme changes in the circulation pattern from the noninteractive 


case. The maximum in rip-current velocity is about 2.6 meters/sec. 
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The failure of the numerical program to compute wave-current inter - 
actions greater than about 50% of the noninteractive case is the inclu- 


sion in the algorithm of a condition which stops the calculation in the 


event that ks0 .- This physically implies that the local wave 


is unable to penetrate beyond this point due to the magnitude of the 
opposing mean current. One possibility to extend computation is to 
assume that when k50, H=0. 


This will obviously yield large 
changes in the spacial derivatives of the H field. The effect of this 
change in the program is presently unknown, but if implemented 
wave diffraction effects should be considered due to possible sharp 
discontinuities in the wave height field. 


3.4 CONCLUSIONS 


1, When the process of wave-current interaction is important, the 
ability for the mean current to convect wave energy across orthogonal 
lines implies that the method of integrating along characteristic lines 
; is no longer valid, Furthermore the wave height can no longer be 

4 simply represented by the product of independent refraction and 


shoaling coefficients, 


2. While only limited interaction results have been obtained thus far, 
it is clear that wave-current interaction is an important physical 
process which can greatly affect incoming wave characteristics with- 
in the nearshore coastal zone. This influence would be especially 
important if a strong circulation pattern exists within the nearshore 


zone. 


3. Since the effects of wave-current interaction on the nearshore 
circulation pattern is so startling, other parallel areas of effort 


i Sle 


should be also considered such as the simultaneous movement of 
bottom material as the current bottom shear stress becomes large, 
A conclusion of a previous study [Noda (1972) | showed that the cir- 
culation pattern and especially the magnitudes of the circulation vel- 
ocities were very sensitive to bottom configurations. The interplay 


of all these factors must be seriously considered. 


4, The preceeding results have opened up interesting areas for fu- 
ture research effort. For the present numerical model some non- 
linearity could be included either from the inertia term or through 
the bottom function term, Also as the wave number becomes k<(0, 


the numerical algorithm should have H- 0, 


Future efforts should be expended ona perturbation expansion tech- 
nique which would allow the use of integrating along characteristic 


lines to obtain approximate solutions for wave-~current interaction, 
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Or on the practical side, these techniques could be directed toward 
an analytic quantification of such basic coastal engineering problems 
as groin spacing and design, 
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4. WAVE AND CURRENT 


4.1 INTRODUCTION 

The question of mass transport or current associated with wave propagation 
in the oceanas well as in channels or rivers has been investigated by nume- 
rous authors, Several types of problems arise in this connection but they 
can be roughly divided into two main classes. The first and most widely 
studied phenomenon is the mass transport due to gravity waves, the second 


is that of wave propagation in the presence of a current. 


The mass transport due to gravity waves has been shown to be of second 

order with respect to wave height. M. L. Dubreil-Jacotin (1934) first proved 
the existence of waves associated with the rotational motion of a perfect 
inviscid fluid and determined that there are an infinity of such solutions 
associated with a more or less arbitrary distribution of vorticity. Since 

the motion is irrotational to the first order, Miche (1944) used the first order 
irrotational salution for a constant finite amplitude two dimensional mcticn 

and calculated the corresponding second order term; the results depend on an 
arbitrary function of depth which can be determined if the vorticity is specified. 
The effect of viscosity was first studied by Longuet-Higgins (1953 and 1960). 
He showed that the boundary layers near the surface and at the bottom were 

of very small extent and that the mass transport velocity above the bottom 
layer and its vertical gradient just below the free surface layer were inde- 
pendent of viscosity. His results coufirmed the existence of vorticity in the 
mean flow and provided boundary conditions which may be used to definc the 
unknown function suggested by Dubreil-Jacotin and used by Miche, Further 
computations of mass transport in cnoidal waves were presented by LeMéhautd 


(1968) who showed that the mass transport is uniform in a vertical plane to the 
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second order of approximation. A further step in the theoretical study of 


mass transport in gravity wave was the study of a spectrum of random waves 


Ming-Shun Chang (1969 and 1970). 


The problem of wave motion in the presence of a current has not been studied 
in the same systematic fashion. However, a number of particular solutions 
to the problem have been presented. For instance, Abdullah (1949) assumes 


an exponentional vorticity distribution, Biesel (1950), a constant vorticity and 


Eliasson and Engelund (1972) a hyperbolic distribution with depth. 


Because of the importance of the latter type of motion, one obvious example 
being the interference of waves and rip current near shore, the assumptions 
and equations used in the last three references will be reviewed, then results 
summarized and an extension of Biegel's solution (1950) to a discrete wave 


spectrum will be presented. 


4.2 BASIC EQUATIONS AND ASSUMPTIONS 
The motion studied here is assumed two dimensional and the fluid is incor- 


pressible. 


The basic coordinate system consists in a horizontal x axis laying along the 
mean water surface line, and a vertical y axis positive upwards. The follow- 


ing notation is used: 


u, Vv x and y components of a particle velocity 
p pressure 
v kinematic viscosity 
; p fluid density 
y stream function 


ua 


gravity constant 
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pasic equations of motion are 


3 u v 

q at aw “ay” Pox oe 
,_ av pg BE Gp Be OP co oe a ew 4.3 

P st * 2 ax * Y ay * B dy 2 eae aa 


The use of these equations imply that any density gradient or Coriolis effects 
have been neglected. Moreover, the terms in parenthesis which will 
eventually be dropped imply a laminar flow, but could be replaced by other 


functions of u or v to approximate turbulent shear stresses. 


Equation (4.1) is identically verified by the introduction of the stream function 


: y such that 
u ay (4. 4) 
vee a (4.5) 


Differentiating (4,2) with respect to y and (4.3) with respect of x and subtract- 
ing and replacing u and v by their values in terms of y yields the funda- 


mental equation for the stream function 


2 2 2 .: 2,2 4, 
or in another familiar form with C= - 4 a y 
d¢ = load . (4.7) 


t » 
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For wave like solutions of wave speed c, y is a function of the variables 


inl semen Go aheeashitanH oI SPIRAEA BER iil on 


{x- tt) and (y}. Assuming a wave spectrum exists with wave numbers 


exte fying from k, to kp, the stream function can be written as 
: 2 
y (x, y, t) -f yo (x - cyt, y) dk (4,8) 


whexs 


c, and y are functions cf k. 
k k 


Equat jon (4, 6) is identically verified by the solution 


v¥ * gy dk = constant. 


The ‘sual irrotational solution is obtained by setting the constant at zero. 


Whew shear stresses are negligible, that is outside some bottom and free 
surfgte boundary layer, the right hand side of (4.6) can be neglected, In this 
regio, the stream function must satisfy the oft used equation: 


2 (9° y) + OY 2 (yg? y) _ avy 2 


2 
= 0, 4, 
at BY 3x See ag eo EO (4.9) 


Asswing that the main effect of shear is to establish initially a non-zero 
vortit ity throughout the fluid but that it can be neglected thcreafter, solutions 
of (4.9) will be sought. The order of magnitude of the neglected term in (4.6) 


will than be evaluated to verify the validity of this approach. 


A sees expansion of the stream function and water height in terms of a 


smal! parameter € will be used. Here ¢€ is a measure of the wave height 
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and it will be assumed that the surface equation 1 (x, t) and se stream 


function ¥ (x, y, t) can be written as 


N 

ie eS (4, 10) 
m = 
N n 
n=9 


where the function % represents the main eifect of the curreht, and will be 
assumed to depend on depth only. Note that here the functions a and Ya 


are of the form (4. 8) that is: 


ue 
1m -f "en FO SE) eR m 21 (4. 8a) 
k 
1 
"e 
t -[ Yan (x -c,t, y) dk n 21 (4, 8b) 
1 


Introducing expression (4.11) for Y into the defining equation (3.9), grouping 
terms in various powers of e¢, aud recalling that by hypothes#s ¥), isa 


function of y only (4.9) becomes 


3 
a¥ ay, ay 
3 2 0 3 2 i 0 an 
e [x (7 ¥,) + ape ise ae ]+ 0") = 0 


ter ¥o (y) (4, 12) 
3 
d¥ oY a-¥ 
re) £3 0 t) 2 1 0 
a 9%) + a Sg) a oo : (4. 13) 
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But irom (4. 8b) 


k 
ae 
Y “f He (> eats yale (4. 14) 
ie 


An additional assumption is now made concerning the functional expression 


of te? that is: 


¥ (x - c,t, y) = <,. (y) cos [k (x - c, t)] ‘ 


Hence, 
5 
¥ -f ¥ (y) cos [k (x - c, t)] dk (4. 14a) 
Ky 
2 “2 [eH ae 
v°¥, = Se ee ae a cos [k(x - c,t)] dk. (4. 14b) 
1 2 k k 
kK) dy 


Introducing (4. 14a) and (4, 14b) in equation (4.13), we get 


This equation will be verified in particular if for all k's inthe interval 


ky xsk< ky 
(= ) ( see ) a ¥, - 
—— -¢ _—_—_ - k y = y * (4. 15) 
dy k dy” k dy” k 
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Hence, given at arbitrary function Yo (y) defining the zeroth order stream 


function due to qa depth dependent current, the components of the first order 


associated wavé stream function are given by equation (4. 15). 


The assumptiony leading to the equation are summarized below: 


{1) There is no density gradient in the fluid. 

{2} Coyiolis forces are neglected. 

(3) Shear stresses are neglected. The validity of this assumption 
wil] be tested for each case considered. 

(4) Wave like solutions only are considered. 

(5) The principal particle velocity component is due to the current 

; arnQ depends on depth only. 

(6) The wave height is small compared to the characteristic length 


if defjning the depth variation of the current. 
4.3 BOUNDARY CONDITIONS 


Note that by definition (equations (4.4) and (4.5)) the stream function is known to 


within an arbitrary function of time. Hence, without loss of generality this 


arbitrary functiun can be set to zero. 


l. The bottom is a line of current. 


Assuming the bottom is defined by 


y = —hb (x) 
oy - cb oF 
- ae (x, t, -h) = “dx ay (x, t, -h) ° 


For h constant this becomes 
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¥ (x, t, -h) = 4 
where % is a constant. 


Replacing ¥ by (4.11) yields the condition 


- = 4, 6 
Lin (-h) ie (4, 16) 
¥. (x, t, -h) = 0 ne i e (4. 17) 
2. The free surface is a line of currer* 


Since the free surface is defined by 
y = nlx, t) 


ay _ an 4 an 
ee ae ie ee y to. 


Replacing y and ¥ by(4.10) and (4.11) and grouping terms in power of e. 
[s* on, an, d¥, 2 
€ x «& t 0) + Be ft) + gS (x, th a 10) +0 (€°) = 0. 


Since Yo is a function of y only, ¥ (0) is a constant. Since, moreover, 

¥ is defined within an arbitrary constant, it is always possible to choose Li (0) = 0 
in which case the value of ¥b is defined by other boundary conditions. Note 

that alternately ty could have been taken as zero and ¥o (0) = + determined 


subsequently. 


The free surface condition hence becomes to first order: 


Yo (0) = 0. (4.18) 
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oY, én, an, at 


eenammaacion + Po aged = = ‘a: 4, 
Sa eet SS Gs O at y 0 (4. 19) 
3. The pressure is known at the surface and is usually assumed 
constant 
From equation (4, 2) and (4.7) 
2 2 2 
ay . ay ate ay aty 1 ap | 
Stay * By Oxdy ~ 3x aye ey. Sy 
and 
vy ay ay ay ay 1 ap 
—— + — = « 
at dx oy 2 ox oy p ay B- 
ox 
Hence, 
a*y an _o’y | ay | oy an avy 
atdy ox dtdx Sy | dxdy ~ 3x ax? 
ay faty an ay |, | an 
ax | 4 2 Ox “Oxdy 8 3x 
y 
1 | op , 3n Op] - 
oa E * 3x oy, ie 
Replacing 1 and ¥ by (4.10) and (4.11), at y=0: 
ey, oy, ay, ay, ay an 
eae A oe a a + 0(e%) = 0 
abby © By Sxdy " Sx 42 a) 
if p(x, 7) = Po a constant. 
Introducing expression (4.19) 
afy dy [ay ay, an an, a*y an 
ati ee es ee ape Od rope | De ge SO at es 
ot oy dy | Oxdy aye Ox ot ay” ox 
(4,20) 
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’ - Now from {4. Ga) 


k 


| 2 
1 -f ui (x -c,t) dk . 
ie \ 


1 


Assuming as for + i {equation (4. 14a) that 


k 
Zz 
1, -[ Mk cos [k(x - c,t)] dk . (4.21) 
ky 


The above boundary conditions will be satisfied if: 
from (4.17) 
¥ (-h) = 0 (4.22) 


from (4. 19) 


[2% |s 
Ye a ee Ty = *O at yy = 0 (4, 23) 
from (4. 20) 
2 
d¥, | ay, d %) qty es 
‘ay fk dy ue dy “k BN = Oat y= 0s 
y 
(4, 24) 


In summary, given a function Yo (y) which describes the effect of a current 

and such that ¥, (0) = 0 the firat order correctionto ¥) describing the 

effect of waves on the stream function is defined asa sy (x, y, t) by equation 
(4.14a). This function must satisfy the differential equation (4.15) with boundary 


conditions described Sy equations (4. 22) through (4, 24). 
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For each function % a change in variabies such that 


i ¥ ok = ¥y - cy (4, 25a) 
i UL = U-cy (4. 25b) 
q XH, 7 x- ot (4. 25e) 
q Vie =v (4, 25e) 
yields the following system of equations: 
3 ar 3 
a 2 Fee | EM fg 
F My See aye ee a Cee) 
dy. Bs dy, 
¥, (-h) = 0. (4. 27) 
fe a ies 
yt dy lk Oat ya Os (4, 28) 
av, oh, av | My 
eer ee ee ME ee + gi, = O aty=0. (4. 29) 
dy, dy dy,” dy. k 


Special solutions of the set of equations (4.15, 4.22, 4.23 and 4, 24) or 
(4. 26 through 4. 29) have been obtained for a single frequency component. 
These solutions are discussed in paragraphs 4. 4 through 4.7 where the 


subscript k has been dropped. In paragraph 4. 8, the extension of these 


cases to a finite wave spectrum will be discussed. 


are 4,4 EXPONENTIAL DISTRIBUTION OF VORTICITY (Abdullah (1949)) 


Take y 
E 2 1 g2 = ay u 
: Sp = V'% =o Wace w u 
° 
corresponding to 
4 ze ay 
q Yo ue 
4 a 
4 = Wt, *Y 
: . ¥y a te: 1) 
4 or from equation (4, 25a) 
E u 
a a ee es ae 
: Yt. a oe (e Ll) -cy. 
4 Note that a is related to % (equation (4. 16}) by 
4 u 
4 a tl 
~ ty ae (e 1). 
The first order solution is given by equation (4, 25) 
ve 2 any et a 2 oe 
=] k° + -—- Ty = 1 k* + oyY. 
da 2 ay u_-c 
y u -c ° 
w 
This equation is Abdullah's equation 11. He gives a solution in terms of a 
power series in ku, for an infinite fluid. The dispersion relation 
for this case is shown to be: 
2 
Ss Ss Ss 
02 02 02 


where 


152 


gir theta ise 


zZ 
kun u u 
. _t w) n _wW n (én + 3) w 
Soi ~ ke{2ntl) . es (ntl) (: )s 3: (nt 1)(2n+3) ( c ) fe | 


2 
ku jn u u 
~ w n(n+1) w n(nt2Z)(4n + 3) ( w 
5o2 ~ Ke(@nth) [a 2! (nt+1) ( ); 3) (nt1)}(2n43) c Ve 


Equation (4. 30) is an implicit relation b tween c and k. It was shown by Abdullah 


(1949), that the effect of the current is to lower the wave propagation speed for 


a fixed frequency. 


u 
Numerical solutions of equation (4.19) for small values of — {smaller than 


. 15) are presented in the cited reference. 


Returning to equation (4.6), the solution given above verifies to order ¢€ that the 


left hand side of the equation is zero. The right hand side is evaluated as 
vvivty) wo [vv(v* ¥)] ~ [yoru]. 


Hence, neglecting the right hand side of equation (4.6) to order ¢ is only consistent 


if varu, is of same order as the terms in 6”, that is = 


2 3 2 * 
E hy as v ¥| ‘ « 7 
2 e 
But since ¥ = ¥ cos k (x - ct) = ¥ cos kX for a single frequency 
2 ay 
2 aa a u_e a 
v*¥ = a¥ - «2 ¥ lcos kx = —“——  ¥ 
1 a 2 ay 
y u_e-7-c¢ 


From (4.28) and (4. 29) 


= dy. 
a ola, 


5 , a¥ = 0 : t 

*s dy ay? 
te oT é 2 
ee oa Qu 
“8 R 2 oO 2 = 2 w 
3 a ‘ (. Yy ax YY )- o [eve hk so, tee! 
ag w 
a +3 
ears 3 q 
ee B20 Bow 2, 
ne 


To relate ¢ to known quantities, note that 


: ay dy 
wp a = — + € wt 
om ba dy dy 
i. < a and the basic assumption is 
; ay ay 
- 1 0 
€ dy << ay ° 
_ In order of magnitude 
eau <<€u ’ 
w w 
Hence, ¢€ is a small parameter for this problem provided ¢ a< < lL. 
| | : ee. seers 
ead From the data presented by Abdullah (1949), for Ms 15 cm/s, a= 10 cm 
. 3 
€<< 10° cm say e€ = 10 cm 
varu = 0 (107% x 107% x 10) = 0 Hor hy sees 


whey §) By comparison 


9 


= 2 r 
0 Meme 10 “x10 k) sec 2 


io) 
R 
< 
- 
iT] 


= 0 io7? k) seer 
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‘ Hence, 0 (v au.) =0 (eau 7k) for wave lengths up to 600 m, and 
4 the assumption that the right hand side of (4.6) can be neglected is justified. 
: 4,5 CONSTANT VORTICITY DISTRIBUTION (Biesel, 1950) 
q Take 
; Cac = ia ¥, = A’ a constant 
a 0 0 : 
. Corresponding to 
5 q uy = Ay + wy, 
j or 7 
a bd 1 - _ + t 
u,. Ay + uy e= Ay + K 
_ A2 
1 sa AM i 


so that ¥y and A are related by 


ey aS 


or calling u, the bottom velocity. 


ws = A (-h) + we 
uc -u 
_ Ww b 

A h e 


weve gd £8 de wheat 


with the boundary conditions (4,27) and (4. 28) 


wins sinh k (y+h) 
¥= - Kn sinh kh : 
Hence, for a single frequency component 


2 , 
¥ - 4 + uy - 6%, K SERRA) cos k (x - ct) +0 ( €*) : 


The velocity components hence become 


| us ay + Ay + ef, Ke SOR TER con k x - ot) +0 (2) 


v= ef Kk SRB EG) gin k(x - ct) + 0 (e%) . 
Introducing the expression for ¥ in (4. 25) 


2 cosh kh A 
. [x seas. 8 | - 8 een 


which is the dispersion relation found by Biesel (1950). 


Since c = 7% where =n is the wave period, Stekes dispersion relation 
ine = g ktanhkh is found for the special case where the vorticity vanishes 
7 (A = 0) and there is no current (uy, = 0). 


For the case of no vorticity but a constant current, the relationship becomes 


w? = g Ktanh kh 
where re is an equivalent wave period defined by 
u 
12 _ ww 
; Ww Ww k ° 


In the general case where the vorticity is given as well as the surface currents, 
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equation {4.31} gives the value of the wave propagation speed as a function 


of wave length. Equation (4.31) can be re-written: 


2 kh ({c-u,) - (c-u,) 
{c-u,) —we 


tanh kh (uy-c) eae 
or 
ee 2 kh (c-u,) ~ (c-u)) ae 
w | tanh kh c-uy et 


Solutions only exist if 


Cd 


(c-u) - (e-u,) | 


(c-uy) [ oe) <gh 


that is 
{c-u.) (c-u,) S gh 


where wy is the bottom (absolute) velocity. 


For a uniform current of speed wu, (that is A = 0) equation (4. 31) yields 


and assuming a realistic wind generated wave where c > a, (Wiegel, 1964) 


and a < uw, 


ec «ec 
u 


so that when there is a vertical velocity gradient the wave speed is lower 
than for a constant current of same surface speed (but not necessarily lower 
than for a constant current at the average current speed), In this case 


0 *(v7¥,) = 0 and 


— 
a“¥ 

v(97¥,) = 9¢ (—+ 
d 
y 


uk 
i=) 


- 7%) cos k (x - ct) 


so that the left hand side of (4.6) is irrelevant for zeroth and first order 


solutions. 
4.6 HYPERBOLIC VORTICITY DISTRIBUTION (Eliasson and Engelund, 1972) 


Here a single frequency component is considered and the following 


relationship is assurned 


C2 0% = vee = BY 


and with boundary conditions (4.27) and (4. 28) 


or 

% = - 4 
Hence, 

Ug = 

Wp (0) = 
or 

K =u 
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é Hence, 

Ys = sinh By 
B 

with 


u, = Ksinh By. 
The first order solution is obtained by solving (4.6) which becomes 


Qa 
SE. ic? + 8%) F 
dy 
with boundary condition (4.27) and (4.28) 


~ i Vic? + 2 | +h 


oa ae 
sinhYk” +B" h 


a | 
rT] 
1 


Hence, 
: Viz 2 
K : ~ 8 (y ) . 
¥= cy+- sinh By - ¢ 9, K Snbvk_+ th cos k (x - ct) + 0 ('e) 
sinhVk* + f* h 
and 


2. pe 
u=ct (uy -¢) cosh By - € "1 KVK7+p% coshVk + +h) cos k (x - ct) 
sinhYk‘“+B* (h) 


+0 (e*) 


: Vie | 2 
v = -et, Kk cer Ei +h sink (x - ct) +0 (¢°) . 
sinhVk" +6" h 


Boundary condition (4.29) gives the following dispersion relation 


peices ge coshVk? +67 h _ 
+B = & (Eliasson and Engelund, 1972). 


sinhVk’ +B” h 


Once again for zeio vorticity (f= 0) and no current, this relation reduces 


2 2 Oe 
to Stokes value where K, B Ber 
2 1 
oo teak, 28 CP ee ea 


For non zero vorticity 


(esatae = amt 5 tankNk’ + p* hb (4, 32) 
Ad Vik +8 

emule Ko +B oh okh 24 

ae = Kk“ +p* h tanh kd 


Hence, the speed of the waves relative to the current is lower than their 


speed in the absence of current. 


In order to evaluate the impact of the viscous terms, consider first the 


expression for ue. The small parameter € is defined by 


| e(u, -e) Vic? + 8? | << |u,| 
or 


6 Vk* + 87 <<«l. 


Then 


vv (v¥5) = 0 (vKB*) = 0 [vu - 8%] 


which needs to be compared to 
[e"¥, y 2 v¥ | = ofe2 K* Vk? +p? K Be] 


= o[e* Vk2 + p2 k Be (a - c)*] 
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‘ Hence, viscous term can be neglected if 


vas 6 V2 +p? k fy -cl 


Taking the dispersion relation into account 


: vBse? +B yeh Ys; ann Kt By : (4. 33) 
4 + 


Pe 
k“ + 8 K* + 6° h 


In most real cases 8 will be relatively small, and for discussing 
equation (4. 33), it will be assumed to be at most of order k. Then (4.33) can 


be approximated by 


vBse - “sgh ae . 


For shallow water kh <1 this becomes 


For example, assuming ¢k = 1072 xa = 107! 
10-7 B 210" 107 h h and po! in centimeters 
a 4 
esyio' nh or py¥k <i0% . 


Assuming § is of the same order as k, v'sco.c iy can be neglected for waves 


length larger than about 5 cm. 


For deep water kh >> 1 


VBS Pos (= 


3 i 
10°? gsio* gh »% afk =< 107? 


so that the same order of magnitude as found above is still valid. 


Hence, viscosity can reasonably be neglected for most of the expected wave 
spectra, It is interesting to compare this case with the case presented in 


paragraph 4,4, 


In both cases the vorticity distribution is exponential in character. The main 
difference is that in the previous case, the solution which was obtained was 
such. that the vorticity was a linear function of the stream function whereas here 
it is a linear function of the stream function relative to axes moving at the wave 
speed. The latter approach was used by Eliasson and Engelund (1972) for a 
single frequency as it yields a closed form solution for the velocity components 
and a fairly simple dispersion relation. However, if a wave spectrum is to be 
analyzed, individual solutions are obtained with respect to axes moving at 


different speeds and the combined solution may be difficult to obtain. 
4.7 GENERAL SOLUTION 


For any given absolute velocity distribution 


ae 


3 


The first order solution is given by 


2a 2 


a“¥ a‘u 

r_ .2 ( i o) ~ 

ae a - ¥- (4, 34) 
2 Yo c dy” 


Boundary conditions (4, 27) and (4, 28) dc “ine the two unknown constants of the 


second order differential equation. 
The dispersion relation is then obtained from (4, 293). 


Equation (4. 34) is singular for Uy=cort © 
and for 


= ~ ik” (u 


97° . 


Since ho | © is not a realistic assumption, the possible motions can be 
obtained either by direct integration of (4,34) as was done in paragraphs 4,5 
and 4.6 or by power series approximation around the two other singular points 
as was done by Abdullah (1949), The dispersion relation in the latter case must 


be solved numerically by succ:ssive approximation, 


It should be noted that in paragraphs 4.4, 4.5 and 4.6 the analyses have assumed 
a single wave length. Because of the non linearity of the basic equations, 
superposition of results cannot be readily performed. Some method of obtain- 
ing solutions for a wave spectrum may then be more valuable to the problem 

of estimating forces due to waves than determining a single component response 


in the presence of complex current distribution. The general discussion of 


paragraphs 4, 2 and 4.3 will provide the means of obtaining this extension. 


4,8 SOLUTIONS FOR A WAVE SPECT"NUM 


l. Consider first the exponencial vorticity distribution. Since % is 


independent of k, as in paragraph 4.4 


u 
yy ee 


and the equation for Y, is from (4. 15) 


ay, 2 au e oy 

See k + ¥ (4 35) 
d é u es c k 

¥ w k 


which is the equation solved by Abdullah (1949). The discussion of 


paragraph 4.4 applies for any wave number k, 


This case corresponds to a wind induced current. Given the wind speed, the 
value uy, can be determined from Wiegel (1964). For each significant 
frequency in the wave spectrum, equations such as (4. 30) can be numerically 
solved for or nomograms could be prepared to obtain c¢, (k). Numerical 
integration of equation (4.35), then (4. 14a) would yield the necessary information 


to compute forces due to the combined action of waves and current. 


2. Constant vorticity 


In this case 


A 2 
¥o = 2 - wy 


Equation (4. 26) gives 


Hence, as in paragraph 4.5 


mie ~ sinh k + h) 
Boe AGS eh Me Spink eh 


and 


Giese ye k cosh kh _ A bs 
w Sk sinh kh ~ ul -c g 
which can readily be solved for, 


Computing the velocity components from ¥ (equation (4, 14a)) and integrating 
over a fixed period, the average velocity components at a given fetch and depth 


are then obiained by a simple numerical quadrative in k. 


3. Hyperbolic vorticity distribution 
For a single frequency it was found that 
¥ = “wi sinh By +c 
where c depended on k according to equation (4.32). Hence, in this case 
Y depends on k, and it is nut possible to use the above expression for ¥ 


0 


where a wave spectrum is considered. 


A possible extension of this solution can be obtained by first defining an 


average velocity ¢ by 


2 
ae rh 1 
o> oe c (k) dk 
2 l k 
1 
and taking 
rae : 
% = B sinhBy + cy. 
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Then from (4. 15) 


~ a¥,, 25 2 x ~ 
[u, ~c) cosh By +c - c,] ae ~ kK’ & | = BY (ul - ¢) cosh By 


As mentioned in paragraph 4.6, the simple equation which is obtained when 


< = c, occurs here at mest for one wave number. In general the equation to be 
solved is 
Fir 2 = 
d he ’ -2 . 8 (uy - c) cosh By ; 
dy" (uy, - @) cosh By tt - cy k 


which is of the form of the equation solved for numerically by Abdullah (1949). 


4. General distribution 

In order to obtain forces and moments on structures placed in a current in the 
presence of waves, it is necessary to compute the velocity components or stream 
function ¥. Fora first order solution, given the wave current effect Xo: the 
stream function ¥ is obtained by evaluating the integral (4, 14a) when ¥, (y) isa 
solution of equation (4,15) and when the functional relationship between c and k 
is usually given by boundary condition (4.24). The process is, therefore, usually 


a lengthy one even with the use of high speed digital computers. 
Results are simplified if ree (y) can be evaluated in closed form. 
From (4,15) it can be seen that 


a™ 
dv¢ aYo/ 3 
ee ie dy 


dy” 20. me ids 
dy k 
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ay 


> + fly) ¥ = 0. 
dy 


Known solutions of such an equation can be obtained for 


a) f(y) = 0 4 = ay tb 
vb) fly) = -¥ ¥ = ae NI engi, SY 
2 -y*/a 
c) fly) = @ntl-y % =e Holy) when n is an integer 
and a is a Hermite Polynomial 
a) fly) = 14 282g 
yee Ny [a Fy) + bI- (y)] 


Jy ly? is a Bessel function and Ais not 


an integer 
or ¥, =v [as (yy) +b YG] if X= a 


is a integer. 


Cases c) and d) already involve functions which must be evaluated numerically 


using a series representation. 


The only fairly simple cases are cases a) and b}. 


For b) 


3 
x2 + 2 ¥o/4¥3 - 42 


a¥, / dy-c) 
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——z = 0 ¥, = ay’ + by since ¥, (0) = 0 


or 


Y= Cy ta fe#Y - . HY] GP aa wie 


Both cases have been studied previously and it was shown that the first one 
only with ¥ independent of k could be extended to a wave spectrum without 


undue complication of the computations. 
4,9 CONCLUSION 


Although the equations required to solve at least up to the first order the 
problem of small amplitude wave propagation in the presence of an arbitrary 
current for an arbitrary wave spectrum have been presented, and a review 
of the special solutions previously obtained for a single wave length has been 


made, the problem in general requires lengthy numerical computations. 


However, for a current whose velocity distribution can be approximated by a 
linear depth dependence, it has been shown that, at most, a single numerical 
quadrature was required to obtain the average velocity components. This 
inethod may then be used to estimate the forces due to wave action in the 
presence of a current. An experimental !nowledge of the current velocity at 


yut a few depths (two minimum) will define the parameters necessary to com- 


pletely solve this problem. 
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APPENDIX A: NOTE ON THE BOTTOM FRICTION APPROXIMATION 
q A dissipative effect imposed by bottom friction comprises an important 
3 term in the momentum equations, The bottom friction is non-linear in 
4 nature, and the coexistence of wave orbital motion and circulation renders 
7 
mathematical formulation of this effect even more complex in the surf 
; zone, 
4 The tangential bottom stress B ofa quadratic form is given by 
4 ey . 4,4 ‘ee 
q Bz=colviv A-l 
q in which cis the friction coefficient, 0 the water density, fe 
3 the resultant velocity vector combining wave orbital and circulatic: 
velocities, e.g. 
+ “ 
Ve Us + (u,v) A-2 
4 > : 
SA in which U, is the wave orbital velocity vector, whose x and y 


components are, respectively; 


4 4 
ve) U_ cos@®@ and-U_ siné6, A-3 
oO oO 


and (u,v) is the circulation velocity vector whose x and y com- 


ponents arc, respectively, u and v, 
» 
P A full expression for v | is then, 
ai 1/2 


Iv | 


[‘w, oad hu es we sin) | 


3 3 2 1/2 
[. +v +2U_u cos@-2U_vsin8+U ] A-4 
re) ° o 


Several methods can be used to simplify Eq. A-4. In the previous inves- 


tigation (Noda, 1972; also Thornton, 1969), two assumptions were made. 
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One was to consider the circulation velocity components u, v_ to be small 


as compared with the wave orbital motion, e.g. 
—-> 
\uU | >> u,v A-5 
fa) 
such that Eq. A-4 reduces to 


lvj= Uy A-6 


Eq. A-l is then rewritten 


4 a 
Beep uU, [U, + tu.v)] As] 


° 
The time average of the bottom stress is 


» n) 
<B>rsep <U, > < [U,+(u,%)] > A-8 


Using the linear theory, the term < u.. > is written as 


. . 2H 
<U, 7 = 3etahkd As 
> 
It was also assumed that the term < [U,, + (u,v)] > may be approximated 
by 
+ + 
<[U,+(u,v)]> = <U,> + < (u,v) > A-10 


~> 
Since < Uy > will vanish over a wave cycle, Eq. A-10 now reduces to 
+ 
<[U, + (u,v)] > = <(u,v) > A-1I] 


Combining Eqs. A-8 and A-11, the time-averaged bottom stress takes the 


form 


t a 


nage DT OR SEIS a 


peg coat 


Noting that the friction terms in the momentum equations are defined by 


> + 
+ <B> <B> 
Foo, = 13 
D(n Fd) pd ae 
then 
2¢U 
° 
F os “u A-14 
x -d 
2c U, 
i = arr: id A-15 


In the present investigation, an attempt was made to carry out a more 
rigorous evaluation of these assumptions, Instead of Eq. A-5, we may 


assume 


> 
JU, |, u >>v A-16 


In other words, we assume that a nearshore circulation contains on- and 
offshore velocity component u which is much longer than the longshore 
component and is not negligible as compared with wave orbital motion, 
This assumption is obviously borne out in a circulation containing a strong 
outflow. Since, because of the refraction, the incidence wave angle in the 


surf zone is generally small, we may further assume 
a@~ 0 A-17 


Using Eqs. A-16 and A-17, Eq. A-4 is simplified as 


4 2 ai l/2 
jv] = [e +a, + 0,7] sie ie A-18 


it? 


pe ree nei ara I ge uenctsd Nlaceaiee, aie 


+ 4 
Using unit vectors i and j inthe x and y _ directions, the bottom 


friction is now rewritten as 


> 
B 


3 > 
ép (u + U,) [u. + (av) | 


ty 
i [ée (u + Us) (u + UL cos 0) 


> 
j [ee (a+ U5) (v- UL sing)] A-19 


+ 


—_) 
The time average of Eq. A-19 is evaluated separately for i and j 


components, e.g. 


> => 4uU 5 
<B> =i+ ép [2 westyd 


> | 2U, 
Le Sp [ae av A-20 


Consequently, the friction terms in the momentum equations are rewritten 


as 
Fo = § [.2 «Ze webu A-21 
FL, § [e+ 2 vu, ov A-Z2 
where 
ae Tsu A-23 


The equation to solve to obtain stream function is of the form [see 
Eq. (2.13) in the text |] 


2 
oy 4 OF oy | 
ae On das eee oy A A-24 
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The choice of the simple friction terms, Eqs, A-14 and A-15, give 


fis f, and f, as functions of x and y only, hence the problem 


reduces to solving a linear differential equation, [see Eq. (2,13) in the 


text], The present computation has used this linear case. 


On the other hand, the choice of the more complex friction terms, Eqs, 

A-21 and A-22, yields fis f, and f, as functions of x,y and ff, 

thus the problem becomes non-linear. Numerical calculation of this 
non-linear case was attempted, but the result exhibited extreme instability, 
hence the lack of convergence. The form of our non-linear equation is more 
complex than any known equations for which the conditions for stable solu- 
tion have been well explored. Further investigation of this problem is 
apparently outside the scope of the present investigation and should be 
reserved for a future study. 


179 


APPENDIX B 
COMPUTER PROGRAMS FOR 


WAVE--CURRENT INTERACTION 


181 


Preceding page blank 


PRAGRAM HATNCINPUT#? 2596, OUTPUT. TAPESS INPUT, FAPEGINUTPLIT, TAPS TESi ee 
#TAPFRSS 12, TAPF9O=SS12P, TAPE IOZ5I2) 


PROGRAM COMPUTES WAVE KINEMATIC AND DYNAMICS BY FINTTE 
DYFFFRENCE INSTEAD OF RY CHARACTERISTICS FRR THE 
GENFRAL CASE OF WAVE eCHPRENT INTERACTING 


AAAs 


4 COMMON UOT 20D VOTO AN I LOT AN PSECTOS PUVA CULT p29 HC TOL20 dy 
q {COCO 19) ,8CH LOD -HRREAK (AY, 19), TACHA III E TOPO I,DONKC TI, AN e 
9 APONYOTH e209) PST (TOe SOI EX( TO, 2D eFYC TO, 20) pwOIN SAO SUN DCT 20)» 
BVELOC O20) LASTOTO,P03,Vi ASTO70~,20) 

a COMMON/CONZ GPL, PT2sRAD EPS OK POV, IX 2,072, Te STRMAGMyN, MM 


4 c 
q c REAN AND WRITE INPUT 
F c 
4 g INTERS0eeeK0) WAVEeCURRENT INTFRACTION, 1) AND V = O,0 
E c TNTE Rep eeewAVERCHRKFNT TNHIFMACTIOM, BFAD U AND V FROM PF CALLED NDAY 
i q C TSTOPVedeeG WIT STORE VELOCITY FIFERS [oem VON PERMANENT FTLES 
: q c TSTORVst-eSTORF U AND V IN PRREMANFIT FITF CALLED NODAV 
: 4 cr TSTNRH=H—=eNnN MOT SINKE 4 ANT THETA FYFL NS UN PERM ANFNY FLEES 
E : is TSTORHS1eeSTORE H AN THETA FIFIDS TN PF CALLED ANDAH 
7 E cr TRE ANHSneeCALCUL ATE THE Ho Anny THETA FIFI eS PTRECTLY 
q E r TREADHS 1 eekF AD THE H OAND THETA FLEIRS Fay PF CALLED BOAR 
: 4 ic TSTARPSheen tT STORE PST FIFLN ONT PRRMAKP NET FILE 
F r TSTWRPSPeeSTRe PST FTE ONTO PRRMANEAT FIER CALTRD SeoINAPST SY 
q Re arpés,1) MeN, TIMAXZ, TIMAXH,TIMAKP, SM, INTE P,TSTORV, IST IRH,TRF ANH, 
i : LTSTORR 
fl F 1 FURMATOCRTYO) 
a E WRITE CH, 10) MeNpTIMAXZ, TTUAKH, TIMAYP YM, TRTER,ISTORV,TSTORM, TREADH 
Le TSTaARP 
10 FOPMAT(IOxX,8TL0/% 
REAN(CS,2) THETA IT, HHA X,OY,T,FPSeA&M,FERICT,PSIMAX 
?P FORMAT CHEIN,?) 
WRITE (He T1) THETAN SHH ,DX,NY,T st PS pate RICT,PSIMAX 
1} FORMATC{OX,LOGte as) 
r 
is COMPUTF ANG DEFINE CONSTANTS 
é 
A=9, ANH? 1 
PTE%.141592P 6536 
PY PrP Tar.) 
RPaMstAN H/rP] 
OxXP=enxe. 0 
NvAsnye2 .O 
STGMAZEPTOsT 
MM toMuM=]{ 
MMTN S44] 
Mysto} 
MP=Moe 
Mpsthed 
Neste 
c CONVERT INT AL VAL TR SOF THE TASNE FR WATER Th RADTANS 
THE TADSTHE TAN RAN 
r abRITF CliNsSTA® FT 
BETTE (HG DEV BoP Te PTA CADIZ Vee STOMA THETAN GL gM Oy 1h? 
LP FARM AT (CTE, ALIAS, TAS) 
Pp 
c INTYTTALTZE VELOCTIVY, THF LA AIR WAVE ETCH TV AMPARLES 
c 
¢ CONVERT DEFER weATER Yate S lif Trth Pat ANI: rth 


AIAN 


50 


55 


51 


3n 
300 


a0 


i a) 


XMAXEDKX AFLOAT (Mel) 
ASYAPT=AMaXMAX 
CALL BVNEM (DS TART, 9,0,0.0,0.0,0,0,RKSTART,A) 
ARE SHSTART ARK START 
TSTARTSASINOSIN(O THE TANI &TANH(CARG)) 
TSTARTSP TS YSTAPRI 
ARG PSARGAD LO 
HSHOHAL SSORT CLO CTANHCARG IS CYL OFARGA/SSINH(ARG2))9) 
HSTARTIHHASOR TC CUSCTHETAN) SCOSCTS TART) 
HSYTARTOHST AKT AHSHITAL 

INITYALTZF FIELPS 
nn sa y2i,me 
7ZSPT tC TSTART@PT)aAFLIAT OTe) SFLOAYOMY) 
HHH= HSTARTRFLOATCTe#1)/FLOAT(MIE) 
Kz NMX#FLNATCT#1) 
SINZSSIN(77} 
cus72.9s(77) 
AH SN J=t,Ne 
UCT, JI= 0.9 
Vit, JI=0,.0 
H(T, J}SHHH 
70T,03222 
SIlTeJIV=ASINZ 
cocre Jv eCosZ 
CONTINUE 
Ao $5 T=2,M 
ni SS J2t, 
CALL OFPTHC(I,J) 
CONTINUE 
no Si T=xMMyN,MM 
ONS AMRAXKFLOAT CT o1) 
NN Sf Jat,n 
N(1T, 7) 200 
CONTYNUF 
REANCIO) COPSICIT, J) ¢Je1eNe),1=1eMM) 
60 TM 400 
DN 3A Y=ei,MM 
PST(T,1)= PSICT,N) 
PSICT,NA)= PST(Y,3) 
CONTINUE 
COANTINGE 
CALL SPEFNR 

FND INITIALIZING “ TFLNS 


READ ' ako V FROM PFRMANENT FILFS IF INTERS1 


TF CINTER (FO, 0) GN TO BO 

RFANCT™) COUCT AI) ,IztieN2),I21¢4M) 
RFANC7) CiVOTe J) 2 JS2eNP D> TH1,MM) 
REWIND 7 

CONTINUE 

Gl 1H APS 

WRITE (6,14) COUCT J) -J=9,10),T=1,) 
FORMAT(//,(10F13,5)) 

WRITECH,1S) CUNCTe J), JEtt MP), Tat ,h) 
FORMAIC//, (9F13.5)) 

WRITE [6,14) COVOT, FT) eJ51-10)-/1=10M) 
WRITF(CH,15) COVOTAI) -JEiigned,Yat,™) 
WRITECA,T4) CCRKCT, J), JE1,19)-T=ieM) 
WRITF(6,15) CORCT, JI -JZLI NA), Fel) 


183 


; ‘ WRITECO,1G) COZCTe Ie JZte AO e TSR? 
F WRITE(H,15) COZCT AD) J2Lt N20 eT ste) 


4 WRITECA,EGYCCSTC Ted) oJtbe iO) e Ted.) 
a: WRITE CO,TSICOST OTS FT) -S2iiehed, tele} 
q WRITE (H,1UVCCCOCT I), dete de Ete 4) 


WRITFCO,ISICCCOCT, J) sJELisN2d,tsie) 
2S COMTINIG 
WRITE (6215) XMAX,DSTART, RESTART, ARC, TS TART, HST ACT 
1% FORMAT (INX, RG4,5) 
; KMAXES 
4 MA On Kol, KMAX 
, BRITE (CA,9T) & 
. 4 OT FAPwWATCIOX,aFRED RACK (YEP FSSA, 77? 
4 Mii Je Yri,e 
4 Mio 7a Jstyne 
§ Tr €® CET. 97> 659 TA 77 
E mE NCL J)2 uct, dd 
ver oc?t.d)s vVet,J) 
: mi TO 7B 
g T7 CORT PNIIE 
MOL OCP, IIS ULaASt(Ts JY 
VPENCP, IIs VEAST CIS) 
TR COMTI 
MHEASTCYT, FT) = HOT, I) 
VEASTOCT, I= VETZI) 
UCP, IIS OLASe Cup Hit, drei aAST(T,I)4 
V(T,II= OLPSe (VUE DCTs J eVELaST (ted) ) 
Te CUMTINE 
TEE AGIS=4 
iz TR CK AER, 173 GE IG 4? 
q RPSie 9,498 
Mn Bt yat,um 
AN $% Jet,Ne 
4 TECARSCLUNOL ACT, J <Ul ASTCT,I)) GTA CEPSIMANSCHLASTOT, J) 
Pike ; TECARSCVIN DCT, dy avi ASTOT, JI) GTC EPSIAARS(VE ASTCT, I) 
7 Be COMTINNE 
TE €TFLAG FA, 09 GH TN BS 
wRTTF (A, 34) FPSU 
$4 FURPMATCINX, «ABOVE TIFRATTON N U AND V CHINVEFRGED FAR ERRORS, 
1F 7,40) 
CALL ¢ xtt 
R7 CONT PME 
BS COUNTING 
WETTFE(CHe14) COUCT,t),fet,10},%st,) 
WRETECH,1S) CCHUOECT, Us Pelt elev, Taste) 
KRTTE (6,19) COVOTe Te Sete 19d eb 1 eM) 
WEITEC HIS) COVE Te J etl i et Abe TSte™) 
Fe CTRF ANH CFG, 9) 6M Th RAS 
EEANCRY CCHET, JS), Jete 12), THe 4) 
REAN(RY CO7CT, J), 351 eNA), PEt pM) 
RFWINP A 
wht TF(6,87) 
AT FURMATOIs//Z, VOX, eVALUIES OF tt AND THETA KEAR FROM PERMANENT FILE Aa, //) 
a : Nir RA Tsi,e 
Nit BR Jat, 
AKG7ZS 7OT,9) 
Sreye,d)e StTATARGZ) 
COT, Ps COS CAPGZ) 
KAR CHIT DSi te 
x j at TN Re 


FLaued 


YE 
VIITEL ALSO 
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i Ia rege 


1 eo ad ie 


PE ree eee 


aon 


ano 


in ee) 


aatino 


85 CONTINUE 
NOW SOLVE FOR THETA 


FALL ANGLE CTTRAX7} 
Patil HE TGHYTCTIMAMH) 
Bo PUNTTNIE 
CAtL SNELL CTHE TAG, HH, AM) 


WRITE OUTPUT 


Gn TH 79 
WRITE CO,60) CONCT se J) eJSte iO) se TaieM) 
60 FORMATCLOX,#WATER DEPTHS ARE w/, C10F13.5)) 
WRITF(H,69) CLOCT ZF) ,JaitsN23,Tsts™) 
61 FURMA.C///,(9F 13.53) 
79 CANTINGE 
WRITE(6262) CCHOT J), 3512102 -1=19) 
G2 FOARMATO//1 0X kVALIES (IF THE WAVE HEIGHT ARF/, (C1GF13.5)) 
WRITECH,61) CCHOT, JD -J2it eed, T3i,) 
WRITECH265) CCIBCI- J), Jet e'ed, Taig) 
63 FORMATOC//10X%,eVALUFS OF THF RREAKING INDEX IRe= ¥8s1 Nf) ARFAKING A 
IND TR=0 RREAKING#/, (19T6)) 
WRYTECH,64) CCHHREAK(C Ts J) ¢J=1-19},T7T=4,"} 
64 FORMAT(C//10X,x*RRFAKING WAVF HETGHI«/,C10F13.5)) 
WRITECH,61) (CHRREAKC ISS) 2, ITIL eNO eTFLeM) 


WRITE VALUES OF DFPTH,H AND THETA FOR SNELLTS LAW REGION 


60 TO 26 
WRITFEC6,709 CCNCTe SJ} eJzl,10)-TSMMINZ MM) 
70 FOPHMATOC//10X,*WATER DEPTHS FOR THE SNEILTS LAW REGTONAaS, (10F13.25)) 
WRITFCHs719 CONC T eS) -JE1t eNO}, 74MIN, MM) 
TL FORMAT(OC///,(09F13.5)) 
WRITECH,72) (CHCIL JD, Sei et0), TZMMIN, MM) 
JT? FORMATC//10X,2WAVE HFIGHTS FOR THE SNELLTS LAW RFGIONA/, CI0FI3,5)) 
WRITECOH,7TL) COHCT SJ) epJEL1ieN2),TSMMIN, MM) 
WRITECG.73) (C20% eS) eJH1,19),T=MMIN, MM) 
73 FORMATOC//L OX, * THETA VALUES FOR THE SNELLTS LAW KFGTON®/, (10F13,5)) 
WRITECH,7TLY CCZ0TA J) > SE1i eNO), TSMMIN, MM) 
26 CONTINUE 
Tr (YSTORH ,FQ, 9) GN TN 8? 
WRTTELR) COHCIL» J) eJS12N2),151,6m) 
WPTTE(B) COZCY, J), S21 -N2),TE1, 4M) 
REWIND 8 
82 CONTINUE 


INITTALIZE THE PST FIFLN 


IF (K ,GE. 2) GU TN 98 
CaLl PSTINTCPSIMAX) 
98 CONTINUE 


SNLVE FOR PSI BY TTFRATION 


CALL STRFAM(TTMAXP,FRICT) 
TF CTSTORP Fa, 9) 6 Tn 96 
WRITFC9) (CCPSICI-J) + JSL¢,N2), 751MM) 
PFWIND 9 

96 TONTINUE 
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oct ant antnancbiruetonds &treendreaee ti it ds cei 


o 


ee ree ee 


am} 
on 


CALt SPEEN 
TF ¢(TSTORV Ea, 4) GN TA al 


WRITECT) 
WRITE CT) 
REWIND 7 
CONTTNUE 
CONTINUE 
STmp 


CCUCTs JV, 521 ,N2), 151.4") 
COVOTs J), J£8,N2), Tet, 44) 
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AAO 


Aas 


ao 


315 


316 


StAPAUTINE STREAMCTIMAXP,FRICT? 


SUBROUTINE COMPITFS THE STREAMLINES GIVEN THE H AND THETA 
FIELDS RY A GAUSS*SEINFt RELAXATION TECHNIQUE ANTH FOR 
WAVE AND NO WA'FeCUPRENT INTERACT TON 


MIMENSTON ANOXC7TO,2C}SDNNY OTS, 20) 

COMMON UIT O,20) ce VETO2 A e ZOET Oe 20) pe ST CTOs ACV COLIN PO ,H(T0- 20), 
1°GCG1,197, 8041519) sHROFAK (41,19)-1R041,193,0070,20),9DDxX(70,20), 
2ODANYC 7TH, 207 ,PS1070,70) -FX(7TO,20) FY CT0,20) aWO70220} 

CUMMOANSCONS GePT,PT2eRAD, FPS,NX, DY NX2,NY2,TsSIGMA,MsNg MM 

FQUIVALENCFE CONDY(125),06(1,1)),CONDYC1e1),-HRREAK(H,1)3) 


STATEMENT FIINCTION 


UPPSICY,J3= CePxSQeeCL,J) + C1, 09F XC TJ) aOxXs2 0} 8PST(I“1,J) + 
1C1L.0 © FXIT,F)eOKX/P II*ePSTC Lei, J) +¢ OXSH*(L,O/DYSQ = FYCIT,J}/DY2)% 
OPSICLT» Jal) & OXSNARCPO/DYSG © FYCT,JI/NY23 ePST CT 2341) /C0IN 

END STATEMENT FIINCT TIONS 

MM) =MMe} 

NiSnet 

MP=N42 

NxSQepxaee 

NYSHEHYee? 

ROXNY=ENxSaspyYSa 

CONS 0411.0 ¢ RAHKNY) 

FYGHT=1,0/8.0 

SIXYNT=1.0/16.n 

FCONS FRICTRSORT (2,085) /P] 

CALCULATE DNDX ANN DADYeesNATFE THAT K,DKDX ANN OKHY ARE 
NOT STORED OUF YO CORE SPACE PROBLEMS 

NO 315 JFienN2 

NNOXCted} = 0.0 

PNOY(1,J)=2 0,0 

CONTINUE 

oO 316 I22,6M 

MO 3lo Jz2,N1 

IJeJe1 

rp=en(IeJJ) 

SS2ST(T,J) 

Ccacn(i1,J) 

CALL WYNUMCDD-CO,SSeUCTe JI +e VCle J) - AKA) 

ARG2= 2, 0eRKe DD 

STNHe=S SINHCARG2) 

TT22 TANHCARG2) 

PuAxs (CUlLeteJet(T=1,J) 3 /Px2 

NuNys (uct edeLeaults,dJet))/P¥2 

AVAx=s (V(Te1,-J)eVifel,J) V/A x2 

OVAY= (VT JtleViTeJet} Y/fy¥? 

NTNOK= (ZCYe1,-Jo7( Tet, J) /0x2 

NIAys (ZC yo lete7(Tedet) )/Py¥2 

FaCse UCToJ ASS = VOES JAC 

FES UCT, JAC + VETS IRSS ¢€ AL SGRARIT.O & APGASSTMH2A)/RK 

DKN¥= RRA(DTINXRF AC = COANUDX = SSaAVDX |= aeDDONX 67, JJI/SINH2]) /FF 

NKOYS RKE(DIDY*FAC = CCOeNDUDY = SS*DVDY = AePNNYII, JJ) /SINHA)/FE 

FF= (1.0 = ARGA/STT2A)/5 INH 

ONOX CT, J)= FFe(RKANDNX(T, JJ) ¢ DDAAKNX) 

DNDY(T,J)2 FERC RKANDDYOIT, JI) + ONeOKDY) 

CONTINUE 

RO 317 Tee2,¥h 
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icant cai 


i: 


AMAIA 


317 


PNDKCis fis LsPiCT,N} 
ONNY CT, 132 GN vit. w) 
ONNOXOC TANASE ANAK, SS 
ONPYCTsNP d= OND. 3? 
CONTINUE 
FND CALCULATION FO DNODX END DANRY 


PRECDMPHTE OFDX/F, DFDY/F AND w 


Nn 300 Is2,mm} 

HN 399 Jee,nt 

JI ses 

PHAN CT, II) 

OPXENAHAHKCOT, J) 

NNV=ENNAYCT, FS) 

Nh Psnneae 

HIZH(T, J} 

Heshtiake 

PHOWS(H( Tet, Jen C lel, Jy} /Nx2 

PHOVE(HCT, Tet eH li, Jeg yo spyve 

DHNXVECHC Tepe Jeroen (Lat, ced paee rete det eH late det SUSI DX PADY?Y 

NHNMXSCHE Tape JI eSeORHET TION (TOL er JIISEXSA 

VENVYES CNC Tp Jeter NeHCT, Tee], J+ P/F YSN 

PINKE(7OTe1 ep JeZtiete ds) snr? 

NINVSC7(TeJetdaZ Ol, Jovy sDYV2 

NIOKVSCLET Eb Jed deoZ (Teta detrvesCretsdeotiel (Tete def VIS (Nx eRDye) 

ATOXXSEZ7 Tepe Fm POF CT, TI4ZCL oie, SIISPYSA 

NINVYESC7 CLs Jaes er a7 Pe TIe7 (IT, Iei VISNYSA 

SSESIC,7) 

C(C2¢ Ve 1 oJ) 

SSPsSSaxP 

CCASCC HK? 

STNPZSTIN(S Nelle J)) 

COSASCNSCEOR7T(T 2 IY) 

CALL WVREUMEIPD OOo SSeNCT eID eVO Te Sg RKGZA) 

AUGPS A, NaRkeDNn 

STNRPSSINHCARG? ) 

aRGis PreNn 

TISTANH CARI) 

TT2= TANHC AKG?) 

AURXS(UCTeL, JJ oulClet,J»IsNKP 

NUDVS (UCT, Pet eul le det dyspye 

OVOVS(V( Tat, JyeveileleJVIsNXP 

NVAVS(VCT, Tel eaevele Jel VI /NVe 

Farsltl(T,J)*#8§ © VEIT, Jerr 

FRIUCT, JIaCO #¢ VOTeIDeSS + Oe5RAR(E 66 & ARGAO/SINHPISRK 

AKERYS PRA (DTOX RF AC = CF AMINY © SSANVDY = ABNDX/STNHDD/EF 

NKAVE RERTNTNYREAC = COeNIPY & SSaDVDY = AeNPY/STNHPI SFE 
NOTFe=e FYCT,ISDENXsE AMD FYCT, DSpe DY se 

MPreHnpenl 

PKes RKAP LO 

AKNNXS HKANDY & NNeDKY 

PKNAYS RKANDY & POLK Y 

FXCT,JVS (DEN xen ARENKAN EST TP Y/R? 4 CNN aN xe? OHH ANDX)/HD 

FYCTe JV = (DENY eRK A RDKIRYSTTOISRK OP 4 (DN ANHNY Oe? Oe RNNY) ZHN 
NOTRee TMENSCGSC 

Fre A Gald od & ARPHRASS TMD) 

COURTS AP, 0eFR © ORS 

CORPS EN © 0,5 

RS FITCH T aH? 


aan 


390 


3906 
391 


392 


393 


36n 


FS= FrueStIne 

Fo: Frectse 

reoc= t,0 + Cf2 

8582 1,0 + S82 

$xX2 AX(CHRERCl? + CHF 2*eSS82) 

Sy¥= Aa€eOrtassPp + COF ear ey 

TAUS SIXINT&H2RFNRSING 

NSXOX=s Re(eFSahT Ox ¢ ONOKEY, JEeECOC} + 2 OADHD ASX /HE 

ASXPYs BAC eFSeOFDY ¢ ONHYCT,JIRCCO) + 2, O*NHDYASK/HY 

MSYNYs RACFSANTOY ¢ ONDYET,J)2SSS) ¢ 2, O8DHDYRSY/HI 

ASYNKsS AR(ESADT OY #€ DANKCT, JY*SSS) + 2, OePMDKRSY/HI 

MTAUXS AFCFCeDTOX ¢ FS*OHDX/HYY &¢ DNDX LT, JPR TALUS FN 

NHPaAyYs AeCKCaDTIDY @ FSeOHOYsHE) & ONDYCI, IT) TAUSEN 

ONDXXS CONAKCTe tp JI-PNDX (Tete, JVI ONE 

BNDYYs (ONDYCIE,J¢1ePNDY CT ,cJel}P/DY2 

ANNOKYs CONNDYC TEAL, JIePNDYCTet,JII/DX2 

SECOND MROFR BERIVITIVES 

DSXAYXS Re CaFS#enTPXVe? ORF ORD TOXRDID eS INZANTDOX®DNDY CE, J} eSIN2* 
IMTNYRONDXCT, J4CCC#ONDXY) + 0 25@HI MDHD Ys CoF SAD TOK4COCADNDX CT, J) 
2e 2, OKDHNXENSXAY/HT ¢ 2, 02S e COHN XY/HL Me NHDOXADKDT/H2) 

NTAUT VS RaCFORNTNVY<=2 NaF SanTHYar?) + O,SeFCeHERNTDYSDHDY @¢ 0, 125% 
LFS*(HI*®DHNYY4DHDYae*2) ¢ RWDNDY( Ts Jia (COS2aDIDY+SIN2RDHDY/HE) + 
CONAY CY s SIAN TAUY/SFN + TAUS(DNDYY/FN © (CONDY(CY 0) /F NJ ane) 

NSYNKYS Bel (2 ORF CADTOVANTNXsFSahTNxy) © STNPENTOYeONDX(T, J) eSItax 
LONONY Crs ZIAD IP X4+SSS*AOND XY) + O0,2S*AHIADHDXACFSANIDY + SSSADNOY(I,J)) 
2¢ PLOADHDYeDSYDX/HI + PL ORSVACDHDXY/HE © DHDYaDenxsHe) 

ATAUXNS Ra CFCaDTOXXe2 Oe F Sahl OXar Py @ O,SeFCeHIRANTOxeOHOY © 0, 1256" 
LFS* CHI*XDHNXX4DHD Kate) + BADNDX CTs» J) *#(COSARDTDX + DHOX*SINZ/HI) + 
PNNOX CTs JIAO TAUXZEN ¢ TAUR (INDY X/SFN © CONDXCT, JD EN) 882) 

FEFCONSHIRSART(CRK/SINH?) SOND 

W(T, J)= GaCCDOSXDYX ¢ NOTAUYY © ASYDXY = DYAUKKIZDN = DDYK(DSXDX + 
IDTAUY)I/ON2 + ON KeCNSYNY «¢ NTAUXISDAZISF 

CONTINUE 

GO 10 360 

WRITFE(6,390) CCPXCT J), 321-19) 5T=1,MM) 

FORMAYC//10X,*#NFDX/FR/,(10613.5)) 

WRITF(6,399) COFXCTe J), JSLisNeds FH 7 MM) 

FORMATO ///, (9613.5) ) 

WRITF(6,392) CCFYCT,J),J21-19)/1T21,MM) 

FORMATC//10X%, aDENY/F#/,(10613.5)) 

WRITF(6,3591) CCFYCT e+ J) -J2iis,Ned,Tai,mMm) 

WRYITE(6,393) (CWC Te Je J1,10}, 131 ,4m) 

FORMATC//10X,2eCT, J) 8@/,0190613,5)) 

WRITE(6,-391) (Ew Te JY, IS11-NA),T=1,4M) 

CONTINUE 


PERFORM ITFRATION FOR PST 


nO 310 TIx1,FTMAXP 

TFLaG=} 

N0 320 [22-41 

PSTNEWZUPPSI(I,2) 

TF CABSCPSINEWePSI(162)).6T. (FRPS#AAS(PSINEW))) TEL AG=O 
PSTCT,2)2PSINEW 

PSTCY,NIJ=PST (Te?) 

PSTNEWRUPPST (I, 3) 

TF CARSCPSINE RePST(1,3)).GCT. CFPS*ARBS(PSTNFW))) TFLAGSO 
PSTCY,3)=PSINEW 

PSTCI,N2=PST(1e33 

MN 850 Jsa,N 


189 


otipbsest tee, 33 


Lm is eas gaan Coats 


PSINEwatiPPSIfi,J) 
TF CAPSUPSYREwWRPSTtT, TI) CGT. FERS aARSERS ENE WI Tht aGes 
PST(TsT=PSINER 

350 COMTINGE 
PST(Tst)SPSirctsn} 

320 TONY TNUE 
JF (TFLAG JER. 42 GE Te sHn 

pO CONTIN 
WRITE(6,330) LIMAXP.FPS,K 

B50 FORMATOS/1OX, @RELAXATION FOR PST FALLED AETERR,I7, 30, aTTERATTOINS 4 
TTTH A REMITPER RRO DF a F it gh, 19K, PKK ES, Fay 
WRETFECA, S359) COPSECT se JV, 329,103, Tot1 Mm) 

$41 FORMATE OY, aL AST VALUFS fF PST AREAS, (ENE TL573 
WRITECA, 352) CLOVSTCT,JD, S29 ,Ne2t},1e1,s) 

BBP FIR MAT OSS /, (C9F 1359) 
nu Tm 39a 

$RA0 wWkttF (6-435) TIT,FRS 

SES FIPMALTE s/s Ox, aREL AXKATTON FOR PST CONVERGED aFTERae, 17, SY,@#TTERAT Libs 
1S aPTR & MAXIMUM ERROR NE aZF IOLA) 
SETTF (A, S36) CCPSTOCI,IY, Foi,igi,T=i,4™)} 

yS5m FORMATOISI OX AT NNVERRED VALUES OF PSI ARE RSs, FTOR 15.573 
wHTYF (6,552) COOSTOL eV, ISILONAV, TEL tM) 

4990 2 Trig 


4 SURRDUTINE ANGLE CITMAKS 


SUBRNUTINE SOLVES FOR THETA BY RELAXATION INCLUD S55 
WAVF «CURRENT INTERACTION 


” 
AMIS 


COMMON /CON/G, PT Pi2sRan, FPS, OX ,DY,DX2o.nY 2 eT SIGMA MN 

COMMON U07TO,20254¥070020),2070, 20) ¢SIETH,203,CO( TH, 20) 2H 70,20), 
1CG0G9,919),8041519) pHBRE AK (41419), TA0G1219),00670,20), 000" (70,20), 
PNOAYCTHsPOVePSICTOS29eFXCTO, AO} SFYC TOe20) 64070220) 


PFRENRM ITERATTON 


OAM 


Nienelt 
N@sNe2 
Mi ame} 
q MosMeg 
E OO POH 1Tz1, TTMAX 
. YFLAG=1 
4 Pn 210 Ils1,4e2 
4 xMelT! 
4 Nii 210 Jz2eNni 
3 CALL NEWANG(T,J,2FL AG) 
an ol 240 COMTINUF 
; TE (TFL AG WFO, 1) GH TN 25 
q 200 CUNTINUE 
a WRITF{6,220) TTHAX 
‘ P20 FORMATCIOX, S3HRELAXATION FOR THETA FAILED AFTER, Th, 3X, 10HTTERATION 
a 1S//) 
“y WRITE(H,221) COZCL ed} e JEL e100, tae) 
221 FORMATCIOx,24HLAST VALUES UF THETA ARF//-(10613.5)) 
WRITF(GH,227) COZCL Fd edeti,s2),7T=1,4) 
22 FIRMAT(///,(9F15.25)) 
CALL FXITT 
P59 WRITF(6e251) ITsEPS 
251 FURMATCIiNY,//7e¢10X%, 3S5HSNLUTION FOR THETA OBTAINED AFTER, T6,3X,435RTT 
1FRATTONS WITH A MAXIMUM RELATIVE FRROR OF,3X,F10,5) 
WRITF(H,2°52) (CZCI eID, JtisiGdstTaiem) 
252 FORMAT(1CX,23HSULUTIONS FUR THETA ARF //,(10F13.5)) 
WRTTFE (4,253) CCOZOT J), SS11,NA) 721, 
253 FURMAT(///,(9F13,5)) 


WRITF THETA IN MFGRFES 


co me ae | 


NO 260 YsiemM 
nN 260 J=i,Ne 
Z7OV, SVZZOT, J eRAN 
P60 CONTINUE 
WRITE(6,2°51) ITs FPS 
WRITE (6,2°S2) CCOZC1 e312, 359,103, TS1,) 
WEITF(6,253) (C2Z0} 23) -JEtieNe),TEis) 
NU 270 Te1," 
AN 270 Js1,Ne 
207, VeZCl eI sPAN 
270 CONTINUE 
QF TURA 


19) 


ie ie ia) 


= 
Reproduced. from 
an available copy 


SUBUNUTIRE SFRANGETs-1, TEL AG) 


SUBROUETME COMPUTFS THE PPOATED ANGLE THETA FOR Tue PFI AXAT TON 
TE CHM Tistht 


CLMMOM PLFA, POV EVE TO PAA FEN PHI ST OPM AD oO UC TO, POP, IN A PO), 


(WCB CAL, 197, SOME, 19), ARF AK CGE, IF) IA CAL LO), 0070,96) ,PONK E10, 20), 


toy 


aho 


AMUN ETO ,ePYPPST C70 e289) EYL TO, eM) EY CTO pe Opn 70,20} 
CRUMON ZCURS GpPTeP TPZ, RADNER SOK OY eDEP NYA e Te STOMA M,N 
STATEME UT Fin TTers 


CHE, FISGCAD MET CT ade FP ECO OE aD FD ELOET Feb eT fet yye O,¢2e54 00701 
Pode be 7 Ter, Seve STC Tat, Pye Cpete TD & CeO) ete Z OP, Jed VeCSiiy, 
PTEevpeST El, lat )d) 

RSET SVSO LAS MEST Late Pde S TC page FeSO be fepeeSt et, dep yyen pasa (C7 et 
Pee dbPa2 Cab, TVR CCO CPt HNC OCP tI) & COLOR THb ye 7 Cte Jeb De Cey, 
Ptetryecif—, 14t3¥) 

Pers Cpe JISC Tape, See T at, yrs? 

Peny( P,IPSCU CT, Jt} sliCt, fet sD? 

P¥OXCT, DPoOVe Lege d.avV Tet, TPPsAK? 

PVOVET? JIE ONGE Te JAI av OT, Pe f Vs? 

COTE Te CORP TS SPECIE? BYP XK, GYOY, GUD AND Mepy CAN UF 
APGTIRDT CAL COD ATER ALD AWHAY STOR i) 

FCT, Vytes ry, Peal ST + vlb, FyaS Et T & yarn ali sa € ARMA LSP Mays Fe 

NROVET DI sOe CE Ny Tar vey, Fy oe STP T#OVEY (Ee SVP - AAP E VOT adel 
LSTMeeDry SEE 

NeNV CL, Ss he lh i STRODE CT, DG) & SENTAOVEXET, TI) = ARPPO ACT Jed ds 
LS ptery/t Fr 

PACLT, FID MET, TIASTNT = UCT, TI ALP ONT 


GREER © AE Feb AT Pera 


rusts ct, ft) 

St'ts SSty,0) 

Sgs te} 

CALL RVR SO CTs FTI CUST,STD VC be FeV le J) eRe A) 

ARR PSP Oe aN GT, TY 

Sp epee Peavey 

RRoRECT,? 

TE (FR COTS O09) tf Th an 

VETTE(H,tS9) Te Je20f eS TV,COST SINT UOT eS VeV ET eT), RK, A 
FL De aT Cynx, a Fh TS HHP AATT yVReeMUTPET Tp dette CIST ASEM Ep iip vg hh pAal, 
1OxX,PTS, 713.5) 
cali o&xTT 
FACTSFAC(CT,J) 

Pep c(STt PeCSteF APE /EF VSD 
PEN PEECUISTASUntar iT VER se 
Chess pNP ae OMERS 

ZF wtCO Stab eoy( tet) = STINT aPe Ae (tT, TY & FCT, Tet abbey = SYet, Ts 
hal ree Ae) 

TF CARBS O/ UF AwmZ ET, TVD SATS CR PSR ARSOAR FADD) Tei arsn 
7iT,JIV=S7"'t 

reer, geese Zee. 4) 

Sper, FPESTS (201, 1) 

TE C.F otk. Pe) Gy FA Aan 

an aa | 

7ET P2701, 2) 

ever,e byscret, Ty) 

Sy, Ur ,PpyesSr( r,t) 

mi TH wa 


~ 


40 


404 


Ane 


TF €F NE. 32 GO TM 4Ot 
NENG? 

7FET,N2A98701,5) 
COCT ANZ ECHO Te J) 
SYTCY,NAISSTCY, 1? 

om mm agg 

TF CJ eNF. ND GU TO 40? 
7OT,1) 5205-8) 

CNT, tyeCN(T,N) 
SUcyey}eSICTeN) 

ef FA age 

MIAN+f 

TF €53 NF. ST} GN To gag 
ZET, 20 =701-35) 

CVC pe2eC0(I,9) 
STCUs2 VEST (12I) 

RETURN 

FRA 
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SUMRIMIT TNE aAVRUM CDS CIS Te SL Ty tig Verh pb) 


SoRROUTINE COMPUTES Tok a AVE RUMP Kee RPT] TP CUD ENG wayk 
CUPREAT Ft PPRACT TIN 


Tet) 3 


COMMONS CES RT, Pi FRAN, FPS EY, OLA, OVD, EY STEM A, yh 
FRSK=0, 00} 
ReaPTo/(TaSERT(G#D}) 
Nis 106 TH1,.450 
AZSTGOMA = LIRR KACHGT = VaR ASTRT 
BPTARKA?P 
ART SWK alt 
FUSEXP CARY) 
Fea1,0/F) 
Stone P O/T TEP A) 
SECH ADEE C Hae? 
TISTANHEARL ) 
Fro fARPKe ET «@ 4/7 
ERS CaeCARGASECHD + TT) + 7 ORC IRAOST + VASTNED A? 
Rik? pas RK w F/T pie 
TR CABRSORKME RERKY Gf F, CARSCERSKARKNE AW) JY Gi Th TEe 
Re SRKENE & 
PHM COMP Ps bee 
yvWTTE (he td) Tp he, Te, V 
POE FORMAT OSS/ TAX GOAT TREAT IOS FOP K PATE RO TO) Fue VERE ARTE R, TH, Ake 
$56513,5) 
cali XM ET 
P19 Un oKeKrnF 
ASSICMA @ PANRACOIST © VaR ASTRT 
TE (Re CFT, 0,0) GO Ti 20 
WETTER (Ae 1 $0) D,COST SLT, eV ,hh yA 
TSA EMUMATC LOY, ARK TS NEGAET VP eefQtiTely D,CONST,STNTU,V ere Ans, 
tiax,7614,93 
CAlL EXIT 
YAN PE TIRA, 


Fen 


¥ 
E 
i 
4 
a 
: 
9 
% 
a 
‘a 
is 
a 


enti 


a Lt a 


aan 


1900 
1001 


SURKNUTINE GROUPCI.J,0CGNxX,DCGNY,FF) 


SUBROUTINE COMPUTES THE GROUP VELOCITY PARAMETFRS CG, 
NEGPX, PCEPY INCLUDING WAVE «CURRENT INTFRACTION 


COMMON UCT, 203-V070220) 62070920) 451070020) pCi 76,20),4N(70,20), 
T7641, 19925041619) HARE AK (OL ,19) TRC 41,193 -h 670,20) ,PDNX70,29), 
ENOAVO TO, 20} ¢PST(70,20) -FXCTOs2N  eFVCTOe20) 9N(70,270) 

COMMAN /CON/S GePTsPI2,RANsEPSsAK PNY eNX2,N¥2sTsSIGMA,M,N 

STATEMENT FUNCTIONS 

NUPKETSFISCUCTe 1, SPH Tet, SPI SAK? 

NUDYCT , JIZCUCT Jeb elit t, Jetyyspye 

AVO¥CT , JIECVCT4e1 Jo oVGTet, Fd) /0K2 

NYNYCTsAJIZCVC Le Iti avi, Jal) I/Py¥2 

NTOXCT es JIS C2CT+t 99) = ZOTHt JVs ¥2 

DIOVCE Te IIe C207 Tet} eZC1,J=t })/0Y¥2 

NOTF IF CORE 1S SUFFICIFNYT DVDX, AVDY,DUDX AND NUDY CAN HE 
APRTIORT CALCULATED AND ARRAY SYNREND 

FCT, JSUCT,JI®COST @ VETA,JIRSINE + O,SeAn(t.0 +4 ARGP/SSTNNZI/SRK 

NKEAX CT oe JIS RKRCCUCT se J) *SINTOVC Te JI eCOST)ADIDX(Y,7) « 

1 (CCOSTXNUDKX CT.) 4 SINTANVNXCIZI)) = ARDDDXCIT, J~4)/SINHZ) /EE 

PKNYCTsJ2& RKACCHCT SI) ASINT = VOT, J aACUST)RDINV(Y,J) «© 
1 CCOSTaABUDVET, J) + SINTENVAYCI,I)I) |= ARDNDY(T, Jat D/SINHO) SFE 

END OF STATEMENT FUNCTIONS 


JJzJ~1 
DF P=n (7, JJ) 
COSTSCN(T, J) 
SINTSST(1,.19) 
CALL WVNUMCDFP,COUSTSSINT UCT» J pVCTe J) ¢RK,A) 
NEXT ORERATTONS COMPUTE THE WAVE RRFAKING HEIGHT 
TASTANH(RK&DEP) 
HBREAK (Te JIZN,12ePTCaTA/SAOK 


COSH1= COSH(RK*DEP) 

SECHSQs 1,0/(CNSH!i*8#2) 

ARG2=2,0*RK*NEP 

SINHP=SSINRC ARG?) 

FNSHASECNSHCARG?) 

SINHSO=SINHP a8? 

FEEt(T.J) 

C=SORT(GRTA/RK) 

FF= 0.58(1.0 + ARGP/SINH?) 

CGOCT-J)= FFeC 

P= CRe(SINH?P |= ARG2SCNSHA) SSINHSY) 

NKNOX=S RKANDOXCL JJ} + NFP KNK OKT, J) 
AKNNYS RKeNNNVCET, JJ) + DEP*eDKNVCT, J) 

A= 0,5*6/(CaRK aw?) 

NCOXZ AXCRKRSECHSARDKODK |= TARNKDOX(1,J)) 
NCNYS ARCRKASFOCHSQADKONY © TARNKDYCI.J)) 
NEOGOX= PeDKOAKX +¢ FFADCNX 

OCGNYs PAaNKNNY + FFaDCOY 

69 ™ 1001 

WRITF(6,1000) T,.J, RK pA,C/DONY,NCDY,DKDNX,DKDNY 
FORMAT(CIOX,219,7615,5) 

CONTINUE 

RF TURN 


éNP 


oe Een I se | 


SURROUTINE DEPTH 1s J) 
SHERDUTINE COMPUTES THE WATER DEPTH AND TTS SPACTAL DERTVITIVES 


COMMON 1O70,2030.V(070,20),7070,20),51070, 20) C070, 703, 4 (75,20), 
ICREAL, PIV ,SCAL, LID HARE AK (COU, IFI, TROOPS I FIGD( 7TH, ANY ~ADAXI70,70), 
EDAPYETOSAOIEPST CTOs PO) PFX CTA GS QO) FFY CTO, 29} HL I0,270) 

COMMON /CON/ GePT es Pl ee RAN, FPS Pe DY MP, NYA, TST CMA, MAN 

ASO, 025 

THIAMAL F/B 0 

Re( On Nee TH PRPs 3,9 

FLAMMaAzaa ln 

AzP0,0 

AL PHAZTSO gf 

ALPHASAL PHASRAN 

YSPXRFINAT( lef} 

YSPVAFLOAT( Je} ) 

TALP AAS TAN CALE HA) 

AKGSEY=X*#TALPRADAPT/SFLAMD A 

SESINC AREY 

Saetaea@ 

Sids89es 

AGE TeX RR THIRD SEs 

PROP XP CARCE ) 

COMIPO SOR AMAARPTAXSFE AND S 

Cernescarsr yj 

ACT, Viocheyel isd &¢ Ate FAS10) 

NAPOV EES JISC aPe ee SOae 

“OM Ed PEA @ ONDVET, FP) e TAL PHA # AMWARER RAST ORCI SO & AREESI,N) 
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Say 


[ie Mae Mae ] 


Maa 


R50 


800 


ROL 


Ao2 


499 


SURROUTING NFWMT CI, J, TFLAG) 


SUAROUTINEG COMPLTFS THE UPDATFO wAVE HEIGHT AND CHECKS FOR 
RRE AK ING 
COMMON UC7O2POV 2 VE7TO22M) o 2070420) e S10 10520) -ENE70, 20) HL 70»26), 
(CO C41 619),5C41,19) -HBREAK (43,19); 18(41,19),0070,20),0DDX(70,20), 
2DDDYL70 420 ¢PSTC70¢20) gFK( T0420) ,FYC 70,20) oWE 70070) 
COMMOK /CON/ GaPTePTesRAN EPS OX DY eDX2,0¥2,7 SIGMA, M,N 


COMPUTE NEW WAVE HFIGHT 


TROT, 523 

NtoNe} 

NP2N4e2 

COLSOVE Ts I) # CCCI, JIESTCI,JVISDY 
COPSCUCT eI) + COCLe J eC 1,5) 70% 
HNFWECCCHAHC Ts Jet) = CO28H(y+1,JIZ00CL = CO2 » §(€T,J}/2,0) 
YF (CHNF eft. HHREAK(CI,J2) GO TH RAO 
HNFW= HRREAK(I, JS} 

IBCI,J}26 

CONTINUE 

TF (ABSCHNEWeH(IT,J)9 .GT. CEPSRARS(CHNENW))) IFLAG=9 
HCT, J) SHNEW 

TF (J .sF. 2) FO TO ROO 
HCIT,NIJSHCT, J} 

TRACIANIJETRCT, J) 

GQ nN agg 

TF (J NE. 3) 60 TO SO! 
HCT»NAPYBHET, J) 

TRCTSNPSTB( Te J) 

GO TO 899 

TF (J .NE. N? GO Th AOP 
W(T,T)2H(T,J) 

TRCT,L)2TBC1,0) 

6 IN 899 

TF (CJ NE. NI) GN TU ADD 

H(T,»? SHOT ,J) 

TR(T,2SIBCI, J) 

RF TURN 

END 
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PLT Ed eee ee 
F re : 


SUBRAUTINE SNELL IE THE TAft, pt, amy 


: 

4 c 

c SURROUTINEG COMPUTES THE SNEULTS Law RAVE HEIGHT, REPTH AND 

q c ORTHAGONA’ ANGULAR NYIRECTION NUISIDE OF THE PERTURIC REACH 

; c 

4 COMMON/TONS CPL, PTA ,RAN FPS DX DY ,UX2.0Y2, Tp STGM Ag MN g MM 

; COMMON 0702 2002 VEFG 201g 7OTO, 200, STC70, 20) CLT O20) HE TN, 20) + 
4 TOGECUL, LVI SCAL LA) eHARFAK (ATID) TECAS SIF APO IN, 209, D0DX(70,20), 
4 PNDOY C70, 20) ¢PSTC70 920) SFXCTO, POI ,E VET 20) pWE THY POY 

a: MMINOMS] 

E NP=ENH2 

; At) 660 PeuMIN, MM 

4 NDsaAMaNx ert oarl Tet) 

4 FALL WVNHMEDN,O,0,0.0,0,9, 0,9, 9K, A) 

! AASRKaDP 

4 ANRSASTINOSTN (THE TAN) ® TANH( AA) ) 

&WNG=P] w Ab‘Q: 


AkG=2,ORAA 
SHOAL ZSSORTCL SOS CTANHC AAD £01 DF ARG/STNACARG)))) 
PEF SESORT(CCOSCTHE TAN) /SOCNSFANG)) 
“4 WHT SHHe SEMA! AREF 
5: SSSSINGCARLY 
rreecns (ata) 
Nr Aho Jt1,he 
OCT,J) 200 
HOT, DY Se VHT 
i 7OT, 3) aL 
E STC1T,J)=55 
F CHOP, JIV=ECC 
HNOAX( Ped )sar 
AMNYCT.,J2=0.8 
600 CONTINUE 
PE TIURM ae 
PNT WD 


vv 
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2 SSR PE EES TRO A TS TREE ey seer eile : ; 
: is) : nn ig Ee gence apg 


SUBROUTINE PSIINTCPSEMAX)} 


he Be J 


SURROUTINEG INTYTTALIZES PSt 


COMMON/SCONS GePI,PT2 cP An, FPS ,DX DY, DK2,DY 2,1 ¢STGMAyg MeN, MM 
COMMON UCTO POV OTOS 20} eZOETONS AGI SIOCTO, 20) ENC TO, 203 eH( 70,280) 
10G6043,19),S8041519) sHRRFAKCG1, 192 IRCOD IFIED TO AO) ,NDPXCING2A0)s 
POPAY ETOAC IsPSTC7TOs29 aE XCTO, 20) FY (70,20),W(70,70) 
NPshe? 
MV {SMMet 
FMsFLOAT(MMY) 
ne 450 JSighe 
PST(C1,J)20,0 
PST(MM,JI20,86 
650 CONTINUE 
AM £60 T=e,MMt 
PSTTSPSIMAX&SINQOPIARLOAT(CT=1)/FM) 
ON 660 Jet,N2 
PSTCT,J)=PSII 
660 CONTINUE 
RETURN 
FAD 


vv 


2 ‘ SUBRNUTING HEIGHT OCTITMAN) 


% c SUBROUTINE COMPHTES THE waVE HEIGHT BY RELAXATION INCLUDING 
a € FFFFCTS OF WAVE=CURRENT INTERACTION 
C 


COAKON UOTO, 20 oVITOs AN s ZETOS ANS STETH s 2GY eCOLTO 62D) sHI7TOe 2035 
ITSTALLIGS, SCO1 419) pHRREAM CGY TFG ER CHT 1F, DOFO, ANY ZPADKE TOAD), 
3 CPHPY CTO, 20s PSTC 70020 pF KCTO SAD oF VOCTOs29} PwOTGS ANY 
q PNEMIIN /CONS GyPT, PT? RAD, EPS, DX, 04,042,072, 1,SIGMA, M,N 


P Mj =Me}{ 
i Nien} 
q ¢ COMPUTE VALUES OF SC€},0) 
3 Mie SOO Ts2.mi 
6G Mt S00 gse,N1 


CALL GROUP CT, J,OCGNK pPLGNY,FF) 
PUBKSCUCT els J eUC Tete JIVs/N¥2 
NUPVYSCUC TL» Teh yetif le Jet ly saya 
MVAXKSEVOET¢1, Deve tate, J Vispye 
OVAYSCVET, Vel devils Tes VISDY? 
PIPXSC7 (Fete Peli ret, J snye 
3 ATAv=S(Z2ZCLe Jel ya7( 1, Jap yvysnve 
Mg SSPa= SICT. J) aar 
e CCAS COT, J) ae 
STRx¥= (2, OebF = O.5}"007 ¢ CRE = N.5)*8SS2 
} SIGYY= (2 ,08FF = 0.573488) + FFF = n,S)atfe 
a TAUXV=S FFAS TCP bP *ACi€ ted) 
. SOT, JI= CHT DP eCSEC Te FVaPIOx | CNC SJ ANTOY) = COUPX + DVPY)  e& 
F. TOCOCT, FaNCGNY & STCT, TAP EDY) © FSTCXXKABUIDY &¢ TANIIXYsDUDY + 
OTAUXYROVAX ¢€ STGYYANVDY) 
60 TA S95 
ARTTE (60590) Je TeCROCT, JIVE GOX, NE BAY ,STOXXs STGYY, TALK Ye SET >} 
SON FuURMATCIXY.215-7615,5) 
SOS COMIT Nite 
S90 CONT] MEF 
‘iPene? 
Me=Me? 


c PERFOR™ TTERATION FOR THE WAVE HEVERT FEL 


Riv GEO TTs1,Me2 
ToMelT 
ry Sao ytst, PrMaY 
TRL AG=I 
Al GaN J=Aenl 
rae: CAlt MFRHTCT > Js TFLAG) 
a4 SPM COMTINtE 
j Te (TFLAG FO, 33 Gt Ti S70 

GRMN COMNTENGIE 

AN TTF (6,540) 1,])7 
540 FORKATCHOX,G9RRELAMATEIO® FOR THE WAVE HETGHT FATI FD Th COVERS, 
1OK,9HWON ROW 12,75 e 5X, SHAFTER, Th, S¥e TP OMTTERATIOAS) 
ERYITE Ch, 509) (HOT) De Jet eed 
Sai FURMATOCNOXY,20FE AST VAllifS OF H APE 4, (10013,5)) 

Cat, FETT 
STA wkIIEF (6,547) 1,ft 
SaP FireeaT Cr OK, asap, AMATI OP FOO WAVE HE TGHMT CONVERGPE O. Rie, Te, dd, 

PSE AFTER, TO, ONT TERATEIINSS) : 
SIO CoN ENE 

RE Tikes 

fen 


“3 *, 


_ 


200 


fe Ber oe Bee 


750 


Ten 
75) 
152 


753 


SUARNUTINE SPEFD 


SUBRAUTIE CUMPEHTFS FHe AND ¥ EPMPONFR IS OF TRE CIRCULATION 
VFLOCTYIFS IN THE PSY FIFLD 


COMMON LLCT OPO SVCTI PAO oe ZITA POV ESTOTN PO) ,COCTIO, 20 ,HECTN 2H), 
TOGCHT,TFI,S CAL, 19D HARE AX CAL, 19 TB C4L ep 19) 2D O76, 20), PONXC70,70), 
AAANYC7O, 20d -PSTC7TOe2O) EY CTO, 2M) pEYCTO,29F WEIN, 20) 
COMMOIN/SCONS GePT,PT2,RAD, EPS, NY PY ,PY eS eDY¥2 2 Ts SIGMA she Ne MM 
MMi SMMe} 

MPN 4+2 

NPONG] 

He FSO Jeeasnt 

uci, 5) 20,4 

VCY- Je 0,0 

tHCMM, JS0.0 

ViMM, JIEePST (MMT, JI COX aN (MM, Jud )) 

AN 750 [22,MM{ 

UCTeI)S -CPSICT, Tet) = PSICT, JL) SONC Tp Jat} epy2) 
V(IT,J)= CPSLCTe1,J) = PSTCTaef ,JDI/SCNC Te Jat aNK2) 
CONTINUE 

NM 760 [=1,MM 

MCT2,2)= UCTLN) 

UCT,N23= HCT,3) 

VCTe1)= VOTLN) 

VCT,N2EVET, 3) 

CONTINUE 

WRITE(6,751) COUCTe J) ,J=Et,tO),IHle6m) 
FORMATO//10X,R11 VELUCITIFS#/.C10F13.5)} 
WRITECO6,752) COUT se JI ¢ JS1i eNO) 2TEts MM) 
FURMATOC///,C9F%3.5)) 

WRITE (6,753) COVCLTeJ),J3J=1,10),T=t 26MM) 
FORMATC//10X,4V VELEUCTTIFS*7,(90F13.5)) 

WRITE (6,752) (OVO eS -JEti NP), T=, MM) 

RETURN 

FANT 


FUNCTION STNM(A) 
F=FXP(A) 

SINHS ce = 1.0/E)/2."9 
RETURA 

FNP 


FUNCTION COSH(A) 
FEYP( aA) 

CNSHE (F + 1,0/6)727,9 
RETURN 

FAN 


201 


1S ef 9,50 5.04 1640 4,an Hn ,on{ U,025 9,94 
gn Pairs) 


ee ee 


i tat 
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AOA 


10 


eo 


PROGRAM MAIN CINPUT SOUTPUT, TAPESZ INPUT, TAPE GEQUTPLT) 


PROGRAM COMPUTES HCX) FROM THE AIERFRENTIAL FOUATION 
FOR ENERGY FOR A WAVE SYSTEM PROPAGATING WITN A CURRENT 
HY RUNGASKUTTA 


PIMENSTON YCL,YP CL) 
COMMON CFL 
RFAN(S,1) NOX,T,D¥ 
FORMATCI10,72F10,2) 
WRITE (6,10) NOX,T,DX 
FORMAT CLOX,#NDK, Ts DXA, 19002F 10,2) 
G=9 RObe21 
PT= 31415976536 
PIAe PT427 5 
CELLS eT /SPTe 

INITIALTZF Y aND YP 
Y(ty= 0.774750 
X= 40,0 
INAxe 0 
CALL RUNES OX,OX%,1,%e¢ YP, INDX) 
WRITF (6,20) Me V¥ CH) SYPCH) 
FORMAT (ION, aXsw, FS ote SK, kHEe,GI565,5%,*0HDX=*, 615.5) 
NM 100 JSt,NNXx 
CALL RUNGS(X,DX,1,Y%e¥P,INDX) 
WRITF(H,7°1) KAYCI D6 YPC1) 
FOPMATOIOX FS ete Mo Gt Sete LOX 2615,5) 
CONT Yntit a 
stop 7 
FN : 
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SURGMUTING RUNGS (Ket Ne ¥eYPRIME, INDE YT) 
NEMENSTON VOL) eYPRIMF CV), ZCt Dewi Pew ASOT de waTy) 


CREPES @ RUNGE KE ETA SIUT TON CF SET WF FIRST ROE PPLE. FORTRAR O99 


A7~AWI7IIA,VAIIAMAAA VMs 


DIMENSTONS Must BF SET FOP FACH PROGRAM 

x TADE PENDENT Val] Ail 

4 TNCPFVEST NEL EA X, MAY GE CHANGER TM VALE 

4) NUMER Rw fF FOUAT PMs 

Y PPPRNDF?T VARTaAnEF at fitk UMP FIMEASTONAL 42H AY 

YRRTME DERIVATIVE HEAT K ORF REMPENS OMAP ARRAY 

Tif PRORKAMMPR MUST SUPPiy THEIPIAL VALERS te ¥Ct? Te YON) 
TAPE X FS A VARTABLE wit} Oe SHEED BE SFT Th 7 Ry RRR FACE 
TRUTT Al Fa TRY Thy Tree Stepriwyter, Joey, The Sel ve @ OIFRERFNE 

SFU OF FQUATTENS OR Ts START SETH FR OINTT PAL CORD TT TONS, 

THE PWOGRAMME MUST wHTTF & SUMPOUTINE CALEER RE RTVE atic Come 
PETES THE NERTVATEVES AND SIMEES Tribe 

THE ARGUMETIT EIST TS SUHIRAUTINE DERIVE (X,tp Y, YER IME } 

TF ¢ TAME X) 5,571 

NO FD Trier 

WI CT) tia YRPREME (CL) 

ZCTIAVOT Ve CuI CTI LS) 

aaYeH/P, 

CALE DFRIVE CAN, Z.YPRIMEY 

NO 5 Teter 

WOCT SH YEOTME (CT) 

ZOTISYO Te eK A(T) 

A=¥eH sa, 

CAtL OFRIVE CAH, 7p YPREME ) 

Mi 4 Paty 

wétT fhe YPRIME CL) 

7OTISZVC Lense} 

A=XeH 

FALL DEWTVE CAN, Z,YPEEMED 

MF Jeger 

waCTJSH*YRRTME (1) 

YOTISVET IGE EC (Sa CwACTIOWSETIVI ENE CT On TIIS&,) 

¥YoaY¥eH 

FALE DERIVE xg, Ye YPRIME) 

nt TP & 

FALL DFRIVE OX ONS Ye YER TIME) 

IAAP X=1 J 
Dp Tink, - 
rep - 
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x : SUBRAUTINE DERIVE OX,%s¥e¥P) 


SURRAUTINE COMPUTES DFRIVITIVES DYDX WHERE YeFCT(x) 


a 


4 AIMENSTON ¥(1),YP(1) 
COMMOe CFL 

Hz 09,9 + OL2ayX 

HuUPy=s AL?2 

q RAD= SNRT(1,0 = 4, OatI/CFL ) 
C2 CFL *C9,5 ¢ O,58R4D) 

: Che 0,5*"6 


# ACMXS eNMlinysRAN 

E POOCNXS O0.5*DC0% 

F YPCHYS YC) amt SeDUDX ¢ NCGNKIZ( A Oe CRC K)) 
q VE TURN 

q END 

4 eno 4.90 -0,10 

4 

q vv 
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